      FUNCTION ABIG ( ARG1, ARG2 )
CD0
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J. D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  FEBRUARY 1987
CD1
CD1   PURPOSE :
CD1
CD1         FINDS THE MAXIMUM VALUE BETWEEN ARG1 AND ARG2
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2     NAME    DEFINITION
CD2
CD2     ARG1     ARGUMENT NUMBER 1
CD2     ARG2     ARGUMENT NUMBER 2
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      ABIG    MAXIMUM VALUE BETWEEN ARG1 AND ARG2
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      NONE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      NONE
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      NONE
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      NONE
C
C***********************************************************************
C
      IF ( ARG1 .GT. ARG2 ) THEN
        ABIG = ARG1
      ELSE
        ABIG = ARG2
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE ADD
CD0
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         ADDS A CONSTRAINT TO THE BASIS. SEE EQNS (13A) AND (13B)
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD2
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD3
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      CHC     PRODUCT OF THE NORMAL OF THE CONSTRAINT AND HDOTC
CD4      CSTARC  PRODUCT OF THE CSTAR OPERATOR AND NORMAL OF THE
CD4              CONSTRAINT
CD4      HDOTC   PRODUCT OF THE H OPERATOR AND NORMAL OF THE CONSTRAINT
CD4      ICNEAR  INDEX OF THE NEAREST CURRENTLY INACTIVE CONSTRAINT
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CSTAR   INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD5              CONSTRAINTS IN THE BASIS. SEE EQN (11)
CD5      H       CAN BE CONSIDERED AS A REDUCED INVERSE HESSIAN OPERATOR
CD5              APPROPRIATE TO THE MANIFOLD FORMED BY THE INTERSECTION
CD5              OF THE CONSTRAINTS. SEE EQN (12)
CD5      IBASIS  INDEX NUMBER O F THE CONSTRAINTS IN THE BASIS
CD5      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD5              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN
CD5              THE BASIS
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      CCDCHC  VALUE OF CK'C/CHC
CD6      HCDCHC  VALUE OF HDOTC/CHC
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONSTROLS AND
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION CSTAR  (MAXX  ,MAXC  )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( CSTAR      ,  CINVRS ( 1,1) )
      EQUIVALENCE ( CHC        ,  COMVMA (   5) )
      EQUIVALENCE ( ICNEAR     ,  KOMVMA (   1) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
      IF (NC .GT. 0) THEN
        DO 20 I=1,NC
          CCDCHC = CSTARC(I)/CHC
          DO 10 J=1,NVAR
            CSTAR(I,J) = CSTAR(I,J) - CCDCHC*HDOTC(J)
   10     CONTINUE
   20   CONTINUE
      ENDIF
C
      NC = NC + 1
      IBASIS(NC) = ICNEAR
C
C...THE NEW CONSTRAINT IS ADD TO THE CSTAR OPERATOR
C
      DO 30 I=1,NVAR
        CSTAR(NC,I) = HDOTC(I)/CHC
   30 CONTINUE
C
C...THE H OPERATOR IS UPDATED
C
      IF (NC .GE. NVAR) THEN
        DO 50 I=1,NVAR
          DO 40 J=1,NVAR
            H(I,J) = ZERO
   40     CONTINUE
   50   CONTINUE
      ELSE
        DO 70 I=1,NVAR
          HCDCHC = HDOTC(I)/CHC
          DO 60 J=1,I
            H(I,J) = H(I,J) - HCDCHC*HDOTC(J)
            H(J,I) = H(I,J)
   60     CONTINUE
   70   CONTINUE
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE BUGAID (IBUG,NUM,XNUM)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  MARCH 1986
CD1
CD1   PURPOSE :
CD1
CD1         PROVIDES EXTRA PRINTOUT FOR CLARIFICATION AND DEBUGGING
CD1         PURPOSES.
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      IBUG    FLAG INDICATING WHICH VARIABLE VALUES ARE TO BE OUTPUT
CD2      NUM     NUMBER OF DATA ARRAY POINTS TO BE PRINTED
CD2      XNUM    VALUE TO BE PRINTED
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      NONE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      NONE
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      IECHO   FILE UNIT NUMBER TO WHICH PRINTOUT IS WRITTEN
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      NONE
C
C***********************************************************************
C
      PARAMETER ( IECHO  =    6 )
C
C***********************************************************************
C
      IF (IBUG .EQ. 1) THEN
C
C...    THE DERIVATIVE OF THE LINE SEARCH FUNCTION IS OUTPUT
C
        WRITE (IECHO,510) XNUM
C
      ELSEIF (IBUG .EQ. 2) THEN
C
C...    THE B MATRIX IS OUTPUT
C
        DO 10 I=1,NUM
          CALL COLUM3 (1,NUM,XDUMMY,I)
   10   CONTINUE
C
      ELSEIF (IBUG .EQ. 3) THEN
C
C...    LINE SEARCH FUNCTION VALUE IS OUTPUT
C
        WRITE (IECHO,530) XNUM
C
      ELSEIF (IBUG .EQ. 4) THEN
C
C...    STEP LENGTH IS OUTPUT
C
        WRITE (IECHO,540) XNUM
C
      ELSEIF (IBUG .EQ. 5) THEN
C
C...    CHANGE IN THE CONTROL VARIABLE (DELTA) IS OUTPUT
C
        CALL COLUM3 (2,NUM,XDUMMY,NDUMMY)
C
      ELSEIF (IBUG .EQ. 6) THEN
C
C...    THE MAXIMUM LAGRANGE MULTIPLIER IS OUTPUT
C
         WRITE (IECHO,560) XNUM
C
       ELSEIF (IBUG .EQ. 7) THEN
C
C...     CURRENT GRADIENT OF THE OBJECTIVE FUNCION (GHAT) IS OUTPUT
C
         CALL COLUM3 (3,NUM,XDUMMY,NDUMMY)
C
      ELSEIF (IBUG .EQ. 8) THEN
C
C...    GRADIENT OF THE LINE SEARCH FUNCTION (GLF) IS OUTPUT
C
        CALL COLUM3 (4,NUM,XDUMMY,NDUMMY)
C
      ELSEIF (IBUG .EQ. 9) THEN
C
C...    VECTOR OF LAGRANGE MULTIPLIERS (VLAM) ARE OUTPUT
C
        CALL COLUM3 (5,NUM,XDUMMY,NDUMMY)
C
      ELSEIF (IBUG .EQ. 10) THEN
C
C...    RMS VALUE OF THE GRADIENT OF THE LINE SEARCH
C...    IS OUTPUT
C
        WRITE (IECHO,570) XNUM
C
      ENDIF
C
C***********************************************************************
C
      RETURN
C
  510 FORMAT (1X,'DFLSA =',1PE16.7)
  530 FORMAT (1X,'FLS =',1PE16.7 )
  540 FORMAT (1X,'STEPL =',1PE16.7 )
  560 FORMAT (1X,'LAMBDA MAX =',1PE16.7 )
  570 FORMAT (1X,'GLFMAG =',1PE16.7 )
C
      END
      SUBROUTINE CHANGP
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         EXCHANGES THE NEW CONSTRAINT WITH THE PASSIVE CONSTRAINT
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      LSIZE   TEMPORARY ARRAY SIZE OF A VECTOR BEING MULTIPLIED
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      AC      PRODUCT OF THE A MATRIX AND CT VECTOR
CD4      CAC     CURVATURE OF THE CONSTRAINT VECTOR CT (= CT*A*C)
CD4      CC      PARTICULAR VALUE OF CSTARC WHICH CORRESPONDS TO
CD4              THE CONSTRAINT BEING REMOVED
CD4      CHC     PRODUCT OF THE NORMAL OF THE CONSTRAINT AND HDOTC
CD4      CKAC    PRODUCT OF THE KTH ROW OF CSTAR ANC AC MATRIX
CD4      CSTAR   INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD4              CONSTRAINTS IN THE BASIS. SEE EQN (11)
CD4      CT      CONSTRAINT VECTOR, CORRESPONDING TO THE LARGEST
CD4              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD4              BASIS
CD4      HDOTC   PRODUCT OF THE H OPERATOR AND NORMAL OF THE CONSTRAINT
CD4      ICNEAR  INDEX OF THE NEAREST CURRENTLY INACTIVE CONSTRAINT
CD4      LAMROW  INDEX OF THE CONSTRAINT, CORRESPONDING TO THE LARGEST
CD4              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD4              BASIS
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN THE
CD4              BASIS
CD4      NVAR    (=N+1) NUMBER OF CONSTROLS PLUS THE ARTIFICIAL VARIABLE
CD4      Y       PRODUCT OF CHC*CAC + CC**2
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CKAC    PRODUCT OF THE NORMAL OF THE CONSTRAINT AND HDOTC
CD5      CSTARC  PRODUCT OF THE CSTAR OPERATOR AND NORMAL OF THE
CD5              CONSTRAINTS
CD5      H       CAN BE CONSIDERED AS A REDUCED INVERSE HESSIAN OPERATOR
CD5              APPROPRIATE TO THE MANIFOLD FORMED BY THE INTERSECTION
CD5              OF THE CONSTRAINTS. SEE EQN (12)
CD5      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS PLUS
CD6              THE ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      UDIVY   VALUE OF U DIVIDED BY Y WHERE U = C*(C'HKC)-HKC(C'C*)
CD6      WDIVY   VALUE W DIVIDED BY Y WHERE  W=HKC(C*'AC*)+C*(C'C*)
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITH CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
C
      DIMENSION CSTAR  (MAXX  ,MAXC  )
      DIMENSION UDIVY  (MAXX  )
      DIMENSION WDIVY  (MAXX  )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( CSTAR      ,  CINVRS ( 1,1) )
      EQUIVALENCE ( CAC        ,  COMVMA (   2) )
      EQUIVALENCE ( CC         ,  COMVMA (   3) )
      EQUIVALENCE ( CHC        ,  COMVMA (   5) )
      EQUIVALENCE ( Y          ,  COMVMA (   9) )
      EQUIVALENCE ( ICNEAR     ,  KOMVMA (   1) )
      EQUIVALENCE ( LAMROW     ,  KOMVMA (   3) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
      LSIZE = 1 + MAXX*(NVAR - 1)
      DO 10 I=1,NC
        CKAC(I) = DOTPRD (NVAR,CSTAR(I,1),LSIZE,MAXX,AC(1),NVAR,1)
   10 CONTINUE
C
      DO 20 I=1,NVAR
        UDIVY(I) = (CHC*CT(I) - CC*HDOTC(I))/Y
        WDIVY(I) = (CAC*HDOTC(I) + CC*CT(I))/Y
   20 CONTINUE
C
C...UPDATE THE H OPERATOR AS PER FORMULA EQN (15B)
C
      DO 40 I=1,NVAR
        DO 30 J=1,I
          H(I,J) = H(I,J) + UDIVY(I)*CT(J) - WDIVY(I)*HDOTC(J)
          H(J,I) = H(I,J)
   30   CONTINUE
   40 CONTINUE
C
C...UPDATE THE CSTAR OPERATOR AS PER FORMULA (15A)
C
      CSTARC (LAMROW) = CSTARC (LAMROW) - FP1
C
      DO 60 I=1,NC
        DO 50 J=1,NVAR
          CSTAR(I,J) = CSTAR(I,J)-CSTARC(I)*WDIVY(J)-CKAC(I)*UDIVY(J)
   50   CONTINUE
   60 CONTINUE
C
      IBASIS(LAMROW) = ICNEAR
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE CHKDEP (DEPEND)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         CHECKS DEPENDENCY OF THE NEW CONSTRAINT
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      DEPEND  FLAG TO INDICATE THE DEPENDENCY OF THE NEW CONSTRAINT
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      AMDAMX  THE MAXIMIMUM LAGRANGE MULTIPLIER
CD4      CAC     CURVATURE OF THE CONSTRAINT VECTOR CT (= CT*A*C)
CD4      CCGCAC  REQUIRED DIRECTION OF SEARCH
CD4      CHC     PRODUCT OF THE NORMAL OF THE CONSTRAINT AND HDOTC
CD4      CSTARC  PRODUCT OF THE CSTAR OPERATOR AND THE NORMAL OF THE
CD4              CONSTRAINT
CD4      GHAT    CURRENT GRADIENT VECTOR OF THE OBJECTIVE FUNCTION.
CD4              SEE EQN (9)
CD4      HDOTC   PRODUCT OF THE H OPERATOR AND THE NORMAL OF THE
CD4              CONSTRAINT
CD4      LAMROW  INDEX OF THE CONSTRAINT, CORRESPONDING TO THE LARGEST
CD4              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD4              BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CC      PARTICULAR VALUE OF CSTARC WHICH CORRESPONDS TO THE
CD5              CONSTRAINT BEING REMOVED
CD5      Y       PRODUCT OF CHC*CAC + CC**2
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      GHC     PRODUCT OF GHAT AND HDOTC
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS PLUS
CD6              THE ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      T1      TEMPORARY PRODUCT TO CHECK DEPENDENCY
CD6      T2      TEMPORARY PRODUCT TO CHECK DEPENDENCY
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM  "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      LOGICAL DEPEND
C
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( AMDAMX     ,  COMVMA (   1) )
      EQUIVALENCE ( CAC        ,  COMVMA (   2) )
      EQUIVALENCE ( CC         ,  COMVMA (   3) )
      EQUIVALENCE ( CCGCAC     ,  COMVMA (   4) )
      EQUIVALENCE ( CHC        ,  COMVMA (   5) )
      EQUIVALENCE ( Y          ,  COMVMA (   9) )
      EQUIVALENCE ( LAMROW     ,  KOMVMA (   3) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
      CC = CSTARC (LAMROW)
      Y = CHC*CAC + CC*CC
      GHC = DOTPRD ( NVAR,GHAT(1),NVAR,1,HDOTC(1),NVAR,1)
      T1 = CCGCAC*Y
      T2 = CHC*(AMDAMX - CCGCAC*CAC) + GHC*CC
C
      IF (T1 .LT. T2) THEN
        DEPEND = .FALSE.
      ENDIF
C
***********************************************************************
C
      RETURN
      END
      SUBROUTINE CKVRTX (NVIOL, SNORM)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         CHECKS THE VERTEX AGAINST THE INACTIVE CONSTRAINTS TO
CD1         MAKE SURE ITS A FEASIBLE POINT
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NVIOL   NUMBER OF VIOLATED CONSTRAINTS
CD3      SNORM   SUM OF THE VIOLATED CONSTRAINT NORMALS
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      BDL     LOWER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1A)
CD4      BDU     UPPER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1A)
CD4      CN      MATRIX OF CONTSTRAINT NORMALS
CD4      D       RIGHT HAND SIDE OF CONSTRAINT EQUATION.
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      IACTIV  STATUS OF EACH OF THE MTOTAL CONSTRAINTS
CD4              =-1 EQUALITY
CD4              = 0 ACTIVE
CD4              = 1 INACTIVE
CD4              = 2 VIOLATED
CD4      MTOTAL  TOTAL NUMBER OF CONSTRAINTS (= M+2*NVAR). INCLUDES THE
CD4              INPUT CONSTRAINTS AND THE UPPER AND LOWER BOUNDARIES
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD4      NVAR2   TWICE THE VALUE OF NVAR. (= 2*(N+1))
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      RESID   PRODUCT OF THE ITH RESIDUAL IN THE BASIS AND RECIPR
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXACT  MAXIMUM SIZE OF ARRAY DEALING WITH IACTIV ARRAY
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6              AND THE ARTIFICIAL VARIABLE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "THE CALCULATION OF FEASIBLE
CD7      POINTS FOR LINEARLY CONSTRAINED OPTIMIZATION PROBLEMS" AERE-
CD7      R.6354 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXACT = MAXX*2 + MAXC )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION SNORM  (MAXC  )
C
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMACT/ RESID  (MAXACT)
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMACT/ IACTIV (MAXACT)
C
      EQUIVALENCE ( MTOTAL     ,  KOMVMA (   4 ) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6 ) )
      EQUIVALENCE ( NVAR2      ,  KOMVMA (   7 ) )
C
C***********************************************************************
C
      NVIOL = 0
      DO 10 I=1,NVAR
        SNORM(I) = ZERO
   10 CONTINUE
C
C...THE CONSTRAINT RESIDUALS, THE NUMBER OF VIOLATED CONSTRAINTS, AND
C...THE SUM OF THEIR NORMALS ARE CALCULATED
C
      DO 30 I=1,MTOTAL
        IF (IACTIV(I) .GE. 1) THEN
C
C...      THE CONSTRAINT IS INACTIVE AND CHECKED FOR VIOLATION
C
          IF (I .LE. NVAR) THEN
            RESID(I) = DELTA(I) - BDL(I)
          ELSEIF (I .LE. NVAR2) THEN
            J = I-NVAR
            RESID(I) = BDU(J) - DELTA(J)
          ELSE
            J = I-NVAR2
            RESID(I) = DOTPRD (NVAR,CN(1,J),NVAR,1,DELTA(1),NVAR,1)
            RESID(I) = RESID(I) -D(J)
          ENDIF
C
          IF (RESID(I) .LT. 0.) THEN
C
C...        THIS CONSTRAINT HAS BEEN VIOLATED
C
            NVIOL = NVIOL + 1
            IACTIV(I) = 2
C
            IF (I .LE. NVAR) THEN
              SNORM(I) = SNORM(I) + FP1
            ELSEIF (I .LE. NVAR2) THEN
              J = I - NVAR
              SNORM(J) = SNORM(J) - FP1
            ELSE
              DO 20 J=1,NVAR
                SNORM(J) = SNORM(J) + CN(J,I-NVAR2)
   20         CONTINUE
            ENDIF
          ENDIF
        ENDIF
   30 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE COLUM3 (IFORM,NUM,XNUM,NUMX)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  MARCH 1986
CD1
CD1   PURPOSE :
CD1
CD1         PROVIDES PRINTOUT DATA IN COLUMNS OF 3.
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      IFORM   FLAG INDICATING WHICH DATA ARRAY TO PRINT
CD2      NUM     NUMBER OF DATA ARRAY POINTS TO BE PRINTED
CD2      NUMX    EXTRA INTEGER NUMBER USED FOR MATRIX COLUMN
CD2              IDENTIFICATION
CD2      XNUM    A DATA ARRAY TO BE PRINTED
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      B       VARIABLE METRIC MATRIX. CORRESPONDS TO THE HESSIAN
CD4              MATRIX IN THE QUADRATIC FUNCTION.
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      GHAT    CURRENT GRADIENT OF THE VECTOR OF THE OBJECTIVE FUNCTION
CD4      GLF     GRADIENT OF THE LAGRANGE FUNCTION. SEE EQNS (3.2)
CD4              AND (3.3)
CD4      VLAM    VECTOR OF LAGRANGE MULTIPLIERS
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      NONE
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FORM    CHARACTER VARIABLE STORING VARIABLE FORMAT STATEMENT
CD6      NFORM   SIZE OF THE FORM ARRAY
CD6      IECHO   FILE UNIT NUMBER TO WHICH PRINTOUT IS WRITTEN
CD6      IEND    LAST POINT OF DATA ARRAY PRINT IN A ROW
CD6      IRPT    NUMBER OF PIECES OF DATA TO PRINT IN A ROW
CD6      ISTART  START POINT OF DATA ARRAY PRINT IN A ROW
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXN    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6      NRPEAT  ARRAY CONTAINING NUMBERS FOR REPETITION FACTORS
CD6      REPEAT  CHARACTER ARRAY CONTAINING REPETITION FACTORS
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      LIMITED TO 72 COLUMNS OF OUTPUT
C
C***********************************************************************
C
      CHARACTER FORM   * 36
      CHARACTER REPEAT *  1
C
      PARAMETER ( IECHO  =    6 )
      PARAMETER ( NFORM  =   10 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXN   = MAXX - 1 )
      PARAMETER ( NRPEAT =    3 )
C
      DIMENSION B      (MAXX  ,MAXX  )
      DIMENSION FORM   (NFORM )
      DIMENSION REPEAT (NRPEAT)
      DIMENSION XNUM   (MAXX  )
C
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAN / GLF    (MAXN  )
     1   ,              GLFA   (MAXN  )
     1   ,              OXA    (MAXN  )
     1   ,              XA     (MAXN  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
C
      EQUIVALENCE ( B          ,  A      ( 1,1) )
C
C***********************************************************************
C
      DATA REPEAT / '1','2','3' /
      DATA FORM   ( 1) /'(Z(1X,2HB(,I2,1H,,I2,2H)=,1PE14.7)) '/
      DATA FORM   ( 2) /'(Z(1X,6HDELTA(,I2,2H)=,1PE13.6))'    /
      DATA FORM   ( 3) /'(Z(1X,5HGHAT(,I2,2H)=,1PE14.7))'     /
      DATA FORM   ( 4) /'(Z(1X,4HGLF(,I2,2H)=,1PE14.7))'      /
      DATA FORM   ( 5) /'(Z(1X,5HVLAM(,I2,2H)=,1PE14.7))'     /
      DATA FORM   ( 6) /'(Z(1X,2HX(,I2,2H)=,1PE16.7))'        /
      DATA FORM   ( 7) /'(Z(1X,2HC(,I2,2H)=,1PE16.7))'        /
      DATA FORM   ( 8) /'(Z(1X,2HG(,I2,2H)=,1PE16.7))'        /
      DATA FORM   ( 9) /'(Z(1X,3HCN(,I2,1H,,I2,2H)=,1PE14.7))'/
C
C***********************************************************************
C
      ISTART = 1
      IF (NUM .LT. 3) THEN
        IRPT = NUM
        IEND = NUM
      ELSE
        IRPT = 3
        IEND = 3
      ENDIF
C
   10 CONTINUE
C
      FORM (IFORM)(2:2) = REPEAT (IRPT)
C
      IF (IFORM .EQ. 1) THEN
        WRITE (IECHO,FORM(1))(I,NUMX,B(I,NUMX),I=ISTART,IEND)
      ELSEIF (IFORM .EQ. 2) THEN
        WRITE (IECHO,FORM(2))(I,DELTA(I),I=ISTART,IEND)
      ELSEIF (IFORM .EQ. 3) THEN
        WRITE (IECHO,FORM(3))(I,GHAT(I),I=ISTART,IEND)
      ELSEIF (IFORM .EQ. 4) THEN
        WRITE (IECHO,FORM(4))(I,GLF(I),I=ISTART,IEND)
      ELSEIF (IFORM .EQ. 5) THEN
        WRITE (IECHO,FORM(5))(I,VLAM(I),I=ISTART,IEND)
      ELSEIF (IFORM .EQ. 9) THEN
        WRITE (IECHO,FORM(9))(I,NUMX,CN(I,NUMX),I=ISTART,IEND)
      ELSE
        WRITE (IECHO,FORM(IFORM))(I,XNUM(I),I=ISTART,IEND)
      ENDIF
C
      IF (IEND .LT. NUM) THEN
        ISTART = IEND + 1
        IF (NUM .LE. IEND+3) THEN
          IEND = NUM
          IRPT = NUM - ISTART + 1
        ELSE
          IEND = IEND + 3
        ENDIF
        GO TO 10
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE DEFER
CD0
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         CALCULATE THE MAXIMUM LAGRANGE MULTIPLIER AND THE GRADIENT
CD1         VECTOR AS THOUGH THE CONSTRAINT IS BEING REMOVED FROM THE
CD1         AUGMENTED BASIS
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD2
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      LSIZE   TEMPORARY ARRAY SIZE FOR A VECTOR BEING MULTIPLIED
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      A       VARIABLE METRIC MATRIX. CORRESPONDS TO THE HESSIAN
CD4              MATRIX IN THE QUADRATIC FUNCTION. SEE EQN (1A)
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      GM      NEGATIVE OF THE GRADIENT VECTOR ( =-G). CORRESPONDS
CD4              TO B IN EQNS (4),(5),(6), AND (8)
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      AMDAMX  THE MAXIMUM LAGRANGE MULTIPLIER
CD5      CT      CONSTRAINT VECTOR, CORRESPONDING TO THE LARGEST
CD5              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FORM THE
CD5              BASIS
CD5      GHAT    CURRENT GRADIENT VECTOR OF THE OBJECTIVE FUNCTION.
CD5              SEE EQN (9)
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS AND
CD6              THE ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( AMDAMX     ,  COMVMA (   1) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
      LSIZE = 1 + MAXX*(NVAR - 1)
      DO 10 I=1,NVAR
        GHAT(I) = DOTPRD(NVAR,A(I,1),LSIZE,MAXX,DELTA(1),NVAR,1) - GM(I)
   10 CONTINUE
C
      AMDAMX = -DOTPRD (NVAR,GHAT(1),NVAR,1,CT(1),NVAR,1)
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE DESIG (MEQ)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         THERE ARE DESIGNATED CONSTRAINTS IN THE BASIS. A FEASIBLE
CD1         VERTEX MUST BE FOUND WITH THE DESIGNATED CONSTRAINTS
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD2
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      KEY     RANK OF THE INVERTED MATRIX
CD3      LSIZE   TEMPORARY ARRAY SIZE FOR A VECTOR BEING MULTIPLIED
CD3      VECDUM  DUMMY VALUE FOR ROUTINE MATIN
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      BDL     LOWER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1A)
CD4      BDU     UPPER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1A)
CD4      CN      MATRIX OF CONSTRAINT NORMALS
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      IBASIS  INDEX NUMBER OF THE CONSTAINTS IN THE BASIS
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN THE
CD4              BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD4      NVAR2   TWICE THE VALUE OF NVAR. (= 2*(N+1))
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CINVRS  INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD5              CONSTRAINTS IN THE BASIS. SEE EQN (2)
CD5      IACTIV  STATUS OF EACH OF THE MTOTAL CONSTRAINTS
CD5              =-1 EQUALITY
CD5              = 0 ACTIVE
CD5              = 1 INACTIVE
CD5              = 2 VIOLATED
CD5      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      CTEMP   TEMPORARY ARRAY TO STORE ROWS OF CINVRS DATA
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      IP      FLAG TO INDENTIFY WHICH ELEMENT OF THE PROJECTION
CD6              IS THE SMALLEST
CD6      IPAFLG  FLAG TO AVOID RECALCULATING P AB-INITIO AT EACH STAGE
CD6              = 0 VECTOR I NOT IN BASIS UPDATE
CD6              = 1 VECTOR I ALREADY IN BASIS UPDATE
CD6      LNEW    INDEX NUMBER OF THE CONSTRAINT BEING BROUGHT INTO
CD6              THE BASIS
CD6      MAXACT  MAXIMUM SIZE OF ARRAY DEALING WITH THE IACTIV ARRAY
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS AND
CD6              THE ARTIFICIAL VARIABLE
CD6      NDESX   COUNTER FOR ADDING REMAINING CONSTRAINTS TO THE
CD6              CINVRS MATRIX
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      P       DIAGONAL ELEMENTS OF THE PROJECTION MATRIX
CD6      PMULT   MULTIPLIER DETERMINED BY WHETHER THE VERTEX IS CLOSER
CD6              TO THE LOWER OR UPPER BOUNDARY
CD6      PSMALL  SMALLEST DIAGONAL ELEMENT OF THE PROJECTION MATRIX
CD6      U       VALUE OF E(I) - V*VPLUS*E(I)
CD6      UUEI    VALUE OF U TRANSPOSE/U TRANSPOSE*E(I)
CD6      V       NORMALS OF THE DESIGNATED CONSTRAINTS
CD6      VPLUSE  VECTOR V+ * E(I)
CD6      VVPLSE  VECTOR V * V+ * E(I)
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "THE CALCULATION OF FEASIBLE
CD7      POINTS FOR LINEARLY CONSTRAINED OPTIMIZATION PROBLEMS" AERE-
CD7      R.6354 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXACT = MAXX*2 + MAXC )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION CTEMP  (MAXX  )
      DIMENSION IPAFLG (MAXX  )
      DIMENSION P      (MAXX  )
      DIMENSION V      (MAXX  ,MAXX  )
      DIMENSION VPLUSE (MAXX  )
      DIMENSION VVPLSE (MAXX  )
C
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMACT/ IACTIV (MAXACT)
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
      EQUIVALENCE ( NVAR2      ,  KOMVMA (   7) )
C
C***********************************************************************
C
C...THE NORMALS V OF THE NC DESIGNATED CONSTRAINTS IN THE BASIS ARE
C...SET UP
C
      DO 30 I=1,NC
        IF (I .LE. MEQ) THEN
          J = IBASIS(I)
          IACTIV(J) = -1
        ELSE
          J = IBASIS(I)
          IACTIV(J) = 0
        ENDIF
C
        IF (IBASIS(I) .LE. NVAR2) THEN
          DO 10 J=1,NVAR
            V(J,I) = ZERO
   10     CONTINUE
C
          IF (IBASIS(I) .LE. NVAR) THEN
C
C...        LOWER BOUNDARY CONSTRAINT IS IN THE BASIS
C
            V(IBASIS(I),I) = FP1
          ELSE
C
C...        UPPER BOUNDARY CONSTRAINT IS IN THE BASIS
C
            JX = IBASIS(I) - NVAR
            V(JX,I) = -FP1
          ENDIF
        ELSE
C
C...      EQUALITY OR INEQUALITY CONSTRAINT OF TYPE EQN (1B) IS IN THE
C...      BASIS
C
          JX = IBASIS(I) - NVAR2
          DO 20 J=1,NVAR
            V(J,I) = CN(J,JX)
   20     CONTINUE
        ENDIF
   30 CONTINUE
C
      IF (NC .EQ. NVAR) THEN
        DO 50 I=1,NVAR
          DO 40 J=1,NVAR
            CINVRS(I,J) = V(I,J)
   40     CONTINUE
   50   CONTINUE
        CALL MATIN(CINVRS,MAXX,NVAR,VECDUM,0,KEY,DETERM)
        IF (KEY .LT. NVAR) THEN
          CALL VMAMSG (15,1,DUMMY)
        ENDIF
        GO TO 999
      ENDIF
C
C...CALCULATE THE INVERSE OF (VTRANPOSE * V) AND STORE IT TEMPORARILY
C...IN CINVRS
C
      DO 60 I=1,NC
        DO 55 J=I,NC
          CINVRS(I,J) = DOTPRD (NVAR,V(1,I),NVAR,1,V(1,J),NVAR,1)
          CINVRS(J,I) = CINVRS (I,J)
   55   CONTINUE
   60 CONTINUE
C
      IF (NC .EQ. 1) THEN
        CINVRS(1,1) = FP1/CINVRS(1,1)
      ELSE
        CALL MATIN (CINVRS,MAXX,NC,VECDUM,0,KEY,DETERM)
        IF (KEY .LT. NC) THEN
          CALL VMAMSG (15,1,DUMMY)
        ENDIF
      ENDIF
C
C...CALCULATE THE GENERALIZED INVERSE OF V (I.E. V+) AND STORE IT IN
C...CINVRS. V+ = (1/(VTRANSPOSE*V)) * VTRANSPOSE
C
      LSIZE = 1 + MAXX*(NC - 1)
      DO 80 I=1,NC
        DO 65 J=1,NC
          CTEMP(J) = CINVRS(I,J)
   65   CONTINUE
        DO 70 J=1,NVAR
          CINVRS(I,J) = DOTPRD(NC,CTEMP,NC,1,V(J,1),LSIZE,MAXX)
   70   CONTINUE
   80 CONTINUE
C
C...CALCULATE THE DIAGONAL ELEMENTS OF THE PROJECTION MATRIX
C...(I.E P(I) = DIAGONAL OF V*V+) AND INITIALIZE BOUND VARIABLES.
C
      DO 90 I=1,NVAR
        P(I) = DOTPRD (NC,V(I,1),LSIZE,MAXX,CINVRS(1,I),NC,1)
        IPAFLG(I) = 0
   90 CONTINUE
      NDESX = NC + 1
C
C...THE REMAINING NVAR-NC CONSTRAINTS ARE ADDED TO THE CINVRS MATRIX
C
      DO 160 LNEW=NDESX,NVAR
        PSMALL=FP1
C
C...    CALCULATE THE SMALLEST BOUND E(I) FROM THE DIAGONAL OF THE
C...    PROJECTION MATRIX AND STORE IT IN PSMALL
C
        DO 100 I=1,NVAR
          IF (IPAFLG(I) .EQ. 0) THEN
            IF (P(I) .LT. PSMALL) THEN
              PSMALL = P(I)
              IP = I
            ENDIF
          ENDIF
  100   CONTINUE
C
        IF (DELTA(IP) - BDL(IP) .GT. BDU(IP) - DELTA(IP)) THEN
          PMULT = -FP1
        ELSE
          PMULT = FP1
        ENDIF
C
C...    CALCULATE THE VECTORS V+ * E(I)
C
        DO 110 I=1,NC
          VPLUSE(I) = PMULT*CINVRS(I,IP)
  110   CONTINUE
C
C...CALCULATE VECTORS V * V+ * E(I)
C
        DO 120 I=1,NVAR
          IF (IPAFLG(I) .EQ. 0) THEN
            VVPLSE(I) = -DOTPRD (NC,V(I,1),LSIZE,MAXX,VPLUSE(1),NC,1)
          ENDIF
          CINVRS(I,IP) = ZERO
  120   CONTINUE
C
C...    CALCULATE U = E(I) - V * V+ * E(I)
C
        U = FP1 + VVPLSE(IP)*PMULT
        IPAFLG(IP) = 1
C
C...    UPDATE V+
C
        DO 140 I=1,NVAR
          IF (IPAFLG(I) .EQ. 0) THEN
            UUEI = VVPLSE(I)/U
            CINVRS (LNEW,I) = UUEI
            DO 130 J=1,NC
              CINVRS(J,I) = CINVRS(J,I) - VPLUSE(J)*UUEI
  130       CONTINUE
          ENDIF
  140   CONTINUE
C
C...    UPDATE THE DIAGONAL ELEMENTS OF THE PROJECTION MATRIX
C
        DO 150 I=1,NVAR
          IF (IPAFLG(I) .EQ. 0) THEN
            P(I) = P(I) + (VVPLSE(I) * VVPLSE(I)) / U
          ENDIF
  150   CONTINUE
C
C...    THE CONSTRAINT BROUGHT INTO THE BASIS IS IDENTIFIED AND INDEXED
C
        CINVRS(LNEW,IP) = PMULT
        IF (PMULT .LT. ZERO) THEN
          IACTIV(IP+NVAR) = 0
          IBASIS(LNEW) = IP + NVAR
        ELSE
          IACTIV(IP) = 0
          IBASIS(LNEW) = IP
        ENDIF
C
        NC = NC + 1
C
  160 CONTINUE
C
  999 CONTINUE
C
C**********************************************************************
C
      RETURN
      END
      FUNCTION DOTPRD (N,VEC1,LSIZV1,INC1,VEC2,LSIZV2,INC2)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1      PERFORMS A DOT PRODUCT OF TWO VECTORS OF LENGTH N.
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      INC1    INCREMENT BETWEEN ELEMENTS OF INPUT VECTOR 1
CD2      INC2    INCREMENT BETWEEN ELEMENTS OF INPUT VECTOR 2
CD2      LSIZV1  TOTAL ARRAY SIZE OF VECTOR 1
CD2      LSIZV2  TOTAL ARRAY SIZE OF VECTOR 2
CD2      N       LENGTH OF THE VECTORS TO BE MULTIPLIED
CD2      VEC1    VECTOR 1
CD2      VEC2    VECTOR 2
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      NONE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      NONE
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      NONE
C
C***********************************************************************
C
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION VEC1   (LSIZV1)
      DIMENSION VEC2   (LSIZV2)
C
C***********************************************************************
C
      DOTPRD = ZERO
C
      DO 10 I=1,N
        DOTPRD = DOTPRD + VEC1 (INC1*(I-1)+1) * VEC2 (INC2*(I-1)+1)
   10 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE EXCHNG
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         UTILIZES THE SIMPLEX METHOD TO EXCHANGE CONSTRAINTS WITHIN
CD1         THE BASIS
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      CSTARC  PRODUCT OF THE CSTAR OPERATOR AND NORMAL OF THE
CD4              CONSTRAINTS
CD4      CT      CONSTRAINT VECTOR, CORRESPONDING TO THE LARGEST
CD4              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD4              BASIS
CD4      ICNEAR  INDEX OF THE NEAREST CURRENTLY INACTIVE CONSTRAINT
CD4      LAMROW  INDEX OF THE CONSTRAINT, CORRESPONDING TO THE LARGEST
CD4              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD4              BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CSTAR   INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD5              CONSTRAINTS IN THE BASIS. SEE EQN (11)
CD5      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS AND
CD6              ARTIFICIAL VARIABLE
CD6      NKOM    KOMVMA ARRAY SIZE
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NKOM   =   15 )
C
      DIMENSION CSTAR  (MAXX  ,MAXC  )
C
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( CSTAR      ,  CINVRS ( 1,1) )
      EQUIVALENCE ( ICNEAR     ,  KOMVMA (   1) )
      EQUIVALENCE ( LAMROW     ,  KOMVMA (   3) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
      RECIP = FP1/CSTARC(LAMROW)
C
      DO 30 I=1,NVAR
        IF (I .EQ. LAMROW) THEN
          DO 10 J=1,NVAR
            CSTAR(I,J) = CSTAR(I,J)*RECIP
   10     CONTINUE
        ELSE
          DO 20 J=1,NVAR
            CSTAR(I,J) = CSTAR(I,J) - RECIP*CSTARC(I)*CT(J)
   20     CONTINUE
        ENDIF
   30 CONTINUE
C
      IBASIS(LAMROW) = ICNEAR
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE GRADLF (G,M,N,IPRINT)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         CALCULATES THE GRADIENT OF THE LAGRANGE FUNCTION
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      G       GRADIENT OF THE OBJECTIVE FUNCTION
CD2      IPRINT  FLAG TO INDICATE AMOUNT OF PRINTOUT
CD2              = 0 NO PRINTOUT
CD2              = 1 VALUES OF X, C, AND F ARE PRINTED
CD2              = 2 DEBUGGING PRINTOUT
CD2      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD2      N       NUMBER OF CONTROL VARIABLES
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      IBUG    FLAG TO INDICATE WHICH DEBUG PRINT IS OUTPUT
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      CN      MATRIX OF CONSTRAINT NORMALS
CD4      VLAM    VECTOR OF LAGRANGE MULTIPLIERS
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      GLF     GRADIENT OF THE LAGRANGE FUNCTION. SEE EQNS (3.2)
CD5              AND (3.3)
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXN    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROLS
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROLS AND THE
CD6              ARTIFICIAL VARIABLE
CD6      GLFMAG  RMS VALUE OF THE GLF VECTOR
CD6      TOL     TOLERANCE FOR EQUALITY ON 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      ALGORITHM AND EQUATIONS CITED BASED ON THE PAPER "A FAST
CD7      ALGORITHM FOR NONLINEARLY CONSTRAINED OPTIMIZATION
CD7      CALCULATIONS" BY M.J.D. POWELL AT THE 1977 DUNDEE
CD7      CONFERENCE ON NUMERICAL ANALYSIS
C
C***********************************************************************
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXN   = MAXX - 1)
      PARAMETER ( TOL    = 1.E-38 )
C
      DIMENSION G      (MAXX  )
C
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAN / GLF    (MAXN  )
     1   ,              GLFA   (MAXN  )
     1   ,              OXA    (MAXN  )
     1   ,              XA     (MAXN  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
C
C***********************************************************************
C
      DO 10 I=1,N
        GLF(I) = G(I)
   10 CONTINUE
C
      IF (M .GE. 1) THEN
        DO 30 I=1,M
          IF (ABS(VLAM(I)) .GT. TOL ) THEN
            DO 20 J=1,N
              GLF(J) =GLF(J) - CN(J,I)*VLAM(I)
   20       CONTINUE
          ENDIF
   30   CONTINUE
      ENDIF
C
      IF (IPRINT .EQ. 2) THEN
        GLFMAG = SQRT ( DOTPRD(N,GLF,MAXN,1,GLF,MAXN,1) )
        IBUG = 10
        CALL BUGAID (IBUG,1,GLFMAG)
        IBUG = 8
        CALL BUGAID (IBUG,N,XDUMMY)
        IF (M .GE. 1) THEN
          IBUG = 9
          CALL BUGAID (IBUG,M,XDUMMY)
        ENDIF
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE INITVM (ISTAT,N,M)
CD0
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         INITIALIZES VARIABLE VALUES NEEDED BY VMACO
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      ISTAT   INDEX TO INDICATE THE STATUS OF THE OPTIMIZATION
CD2              ALGORITHM.
CD2              =-2 INITIAL PASS THROUGH VMACO WITH THE B MATRIX
CD2                  ALREADY INITIALIZED
CD2              =-1 INITIAL PASS THROUGH VMACO
CD2              = 0 CONTINUE ALGORITHM CALCULATIONS
CD2              = 1 CONVERGE WITHIN REQUIRED ACURRACY
CD2              > 1 ERROR ENCOUNTERED
CD2      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD2      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD2      N       NUMBER OF CONTROL VARIABLES
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      IBUGFG  FLAG TO INDICATE EXTRA PRINTOUT FOR DEBUGGING PURPOSES
CD3
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      NONE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      B       VARIABLE METRIC MATRIX. CORRESPONDS TO THE HESSIAN
CD5              MATRIX IN THE QUADRATIC FUNCTION. SEE EQN (2.1)
CD5      ITER    COUNTER FOR NUMBER OF ITERATIONS COMPLETED
CD5      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5      NVMA    NUMBER OF CALLS TO VMACO
CD5      VMU     MULTIPLIERS FOR THE LINE SEARCH FUNCTION. EQN(4.2)
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS PLUS
CD6              THE ARTIFICIAL VARIABLE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7      EQUATIONS CITED FROM "A FAST ALGORITHM FOR NONLINEARLY
CD7      CONSTRAINED OPTIMIZATION CALCULATIONS" BY M.J.D. POWELL
C
C***********************************************************************
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION B      (MAXX  ,MAXX  )
C
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( B          ,  A      ( 1,1) )
      EQUIVALENCE ( ITER       ,  KOMVMA (   2) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
      EQUIVALENCE ( NVMA       ,  KOMVMA (   8) )
      EQUIVALENCE ( IBUGFG     ,  KOMVMA (  11) )
C
C***********************************************************************
C
      IBUGFG = 0
      ITER = 0
      NVMA = 1
      NVAR = N + 1
C
      IF (ISTAT .EQ. -1) THEN
C
C...    INITIALIZE THE B MATRIX
C
        DO 20 I=1,N
          DO 10 J=1,N
            B(I,J)= ZERO
  10      CONTINUE
          B(I,I)= FP1
  20    CONTINUE
      ENDIF
C
      IF (M .GE. 1) THEN
        DO 30 I=1,M
          VMU(I) = ZERO
   30   CONTINUE
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE LGRANG
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         SETS UP THE LAGRANGIAN FUNCTION FOR THE QUADRATIC
CD1         PROGRAMMING ITERATION
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      LSIZE   TEMPORARY ARRAY SIZE OF A VECTOR BEING MULTIPLIED
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      A       VARIABLE METRIC MATRIX. CORRESPONDS TO THE HESSIAN
CD4              MATRIX IN THE QUADRATIC FUNCTION. SEE EQN (1A)
CD4              EQUIVALENT TO THE B MATRIX IN VMACO
CD4      CSTAR   INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD4              CONSTRAINTS IN THE BASIS. SEE EQN (11)
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      GM      NEGATIVE OF THE GRADIENT VECTOR (= -G). CORRESPONDS
CD4              TO B IN EQNS (4),(5),(6), AND (8)
CD4      IACTIV  STATUS OF EACH OF THE MTOTAL CONSTAINTS IN THE BASIS
CD4              =-1 EQUALITY
CD4              = 0 ACTIVE
CD4              = 1 INACTIVE
CD4              = 2 VIOLATED
CD4      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN THE
CD4              BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      AMDAMX  THE MAXIMUM LAGRANGE MULTIPLIER
CD5      GHAT    CURRENT GRADIENT VECTOR OF THE OBJECTIVE FUNCTION.
CD5              SEE EQN (9)
CD5      LAMROW  INDEX OF THE CONSTRAINT, CORRESPONDING TO THE LARGEST
CD5              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM
CD5              THE BASIS
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      AMBDA   LAGRANGE MULTIPLIER (LAMBDA) (= -CSTAR*GHAT)
CD6      MAXACT  MAXIMUM SIZE OF ARRAY DEALING WITH THE IACTIV ARRAY
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROLS AND THE
CD6              ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      XINF    NUMBER FOR INFINITY
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXACT = MAXX*2 + MAXC )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( XINF   = 1.E38 )
C
      DIMENSION CSTAR  (MAXX  ,MAXC  )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMACT/ IACTIV (MAXACT)
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( CSTAR      ,  CINVRS ( 1,1) )
      EQUIVALENCE ( AMDAMX     ,  COMVMA (   1) )
      EQUIVALENCE ( LAMROW     ,  KOMVMA (   3) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
C...THE GRADIENT VECTOR OF THE OBJECTIVE FUNCTION IS CALCULATED
C
      LSIZE = 1 + MAXX*(NVAR-1)
      DO 10 I=1,NVAR
        GHAT(I) = DOTPRD(NVAR,A(I,1),LSIZE,MAXX,DELTA(1),NVAR,1) - GM(I)
   10 CONTINUE
C
C...THE MAXIMUM LAGRANGE MULTIPLIER IS CALCULATED
C
      AMDAMX = -XINF
C
      DO 20 I=1,NC
        IF (IACTIV(IBASIS(I)) .NE. -1) THEN
          AMBDA = -DOTPRD (NVAR,CSTAR(I,1),LSIZE,MAXX,GHAT(1),NVAR,1)
          IF (AMBDA .GT. AMDAMX) THEN
            AMDAMX = AMBDA
            LAMROW = I
          ENDIF
        ENDIF
   20 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE LNSRCH (C,F,G,ISTAT,M,MAXFUN,MEQ,N,NVMITR,NVMLIN,
     1                   X,BFLAG,IPRINT)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         CALCULATES THE NECESSARY INFORMATION FOR THE LINE SEARCH
CD1         OBJECTIVE FUNCTION
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      C       VECTOR OF CONSTRAINT VARIABLES
CD2      F       VALUE OF THE OBJECTIVE FUNCTION (OPTIMIZED VARIABLE)
CD2      G       GRADIENT OF THE OBJECTIVE FUNCTION
CD2      IPRINT  FLAG TO INDICATE AMOUNT OF PRINTOUT
CD2              = 0 NO PRINTOUT
CD2              = 1 VALUES OF X, C, AND F ARE PRINTED
CD2              = 2 DEBUGGING PRINTOUT
CD2      ISTAT   INDEX TO INDICATE THE STATUS OF THE OPTIMIZATION
CD2              VARIABLE
CD2              =-1 INITIAL PASS THROUGH VMACO
CD2              = 0 CONTINUE ALGORITH CALCULATIONS
CD2              = 1 CONVERGE WITHIN REQUIRED ACCURACY
CD2              > 1 ERROR ENCOUNTERED
CD2      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD2      MAXFUN  MAXIMUM NUMBER OF TIMES VMACO CAN BE CALLED
CD2      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD2      N       NUMBER OF CONTROL VARIABLES
CD2      NVMITR  VALUE OF NVMA AT THE START OF AN ITERATION
CD2      NVMLIN  CONTROLS THE ERROR RETURN FROM THE LINE SEARCH OBJEC-
CD2              TIVE FUNCTION. EQN(4.2)
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      BFLAG   FLAG FOR WHETHER THE B MATRIX SHOULD BE UPDATED
CD3      G       GRADIENT OF THE OBJECTIVE FUNCTION
CD3      IBUG    FLAG TO INDICATED WHICH DEGUG OUTPUT IS TO PRINTED
CD3      ISTAT   INDEX TO INDICATE THE STATUS OF THE OPTIMIZATION
CD3              ALGORITHM
CD3              =-2 INITIAL PASS THROUGH VMACO WITH THE B MATRIX
CD3                  ALREADY INITIALIZED
CD3              =-1 INITIAL PASS THROUGH VMACO
CD3              = 0 CONTINUE ALGORITH CALCULATIONS
CD3              = 1 CONVERGE WITHIN REQUIRED ACCURACY
CD3              > 1 ERROR ENCOUNTERED
CD3      MAXFUN  MAXIMUM NUMBER OF TIMES VMACO CAN BE CALLED
CD4      STOFLG  FLAG FOR WHETHER THE ROUTINE STOVMU SHOULD BE CALLED
CD4      X       VECTOR OF CONTROL VARIABLES
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH
CD4      GDELTA  SCALAR PRODUCT OF GRADIENT OF THE OBJECTIVE FUNCTION
CD4              AND DELTA
CD4      NVMA    NUMBER OF CALLS TO VMACO
CD4      VMU     MULTIPLIERS FOR THE LINE SEARCH FUNCTION. EQN(4.2)
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      DFLSAV  SAVED VALUE OF DFLSA
CD5      FLSA    INITIAL VALUE OF FLS. CORRESPONDS TO PHI(0) IN EQN(4.7)
CD5      OXA     SAVED VALUE OF THE VARIABLE XA
CD5      XA      INITIAL GUESS AT CONTROL VARIBLE VALUES
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      CHORD   SLOPE OF THE ARMIJO CHORD
CD6      DFLSA   DERIVATIVE OF FLS. CORRESPONDS TO CAPITAL DELTA OR
CD6              PHI'(0) IN EQN (4.7)
CD6      DIFFLS  DIFFERENCE BETWEEN THE CURRENT LINE SEARCH OBJECTIVE
CD6              FUNCTION AND THE PREVIOUS VALUE. EQN (4.8)
CD6      FLS     VALUE OF THE LINE SEARCH OBJECTIVE FUNCTION. EQN (4.2)
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      HALF    FLOATING POINT NUMBER 0.5
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROLS PLUS THE
CD6              ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      STEPL   STEP-LENGTH FOR THE SEARCH DIRECTION. CORRESPONDS TO
CD6              ALPHA IN EQN (2.2)
CD6      STPMAX  LIMITS THE REDUCTION IN THE LINE SEARCH STEP-LENGTH
CD6      SUMINF  WEIGHTED SUM OF INFEASIBILITIES. EQN (4.2)
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS CITED FROM "A FAST ALGORITHM FOR NONLINEARLY
CD7      CONSTRAINED OPTIMIZATION CALCULATIONS" BY M.J.D. POWELL
C
C***********************************************************************
C
      LOGICAL BFLAG
C
      PARAMETER ( CHORD  =  0.1 )
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( HALF   =  0.5 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXN   = MAXX - 1)
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( STPMAX =  0.1 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION C      (MAXC  )
      DIMENSION G      (MAXX  )
      DIMENSION X      (MAXX  )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAN / GLF    (MAXN  )
     1   ,              GLFA   (MAXN  )
     1   ,              OXA    (MAXN  )
     1   ,              XA     (MAXN  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( GDELTA     ,  COMVMA (   7) )
      EQUIVALENCE ( DFLSAV     ,  COMVMA (  11) )
      EQUIVALENCE ( FLSA       ,  COMVMA (  12) )
      EQUIVALENCE ( NVMA       ,  KOMVMA (   8) )
C
C***********************************************************************
C
      SUMINF = ZERO
C
      IF (M .GT. 0) THEN
        DO 10 I=1,M
          IF (I .LE. MEQ) THEN
            SUMINF = SUMINF + VMU(I)*ABS(C(I))
          ELSE
            SUMINF = SUMINF + VMU(I)*ABIG(ZERO,-C(I))
          ENDIF
   10   CONTINUE
      ENDIF
C
      FLS = F + SUMINF
C
      IF (NVMA .EQ. NVMITR) THEN
C
C...    SET THE INITIAL CONDITIONS FOR THE LINE SEARCH
C
        FLSA = FLS
        DFLSA = GDELTA - DELTA(N+1)*SUMINF
C
        IF (DFLSA .GE. ZERO) THEN
          ISTAT = 7
          CALL VMAMSG(ISTAT,0,ZERO)
          GO TO 999
        ELSE
          STEPL = FP1
C
C...      CALCULATE THE NEW CONTROL VARIABLE GUESSES AND SAVE THE OLD
C
          DO 20 I=1,N
            XA(I) = OXA(I)
            X(I) = XA(I) + DELTA(I)
   20     CONTINUE
C
          DFLSA = STEPL*DFLSA
          DFLSAV = DFLSA
        ENDIF
C
        IF (NVMA .GE. MAXFUN) THEN
          DO 30 I=1,N
            X(I) = XA(I)
   30     CONTINUE
C
          ISTAT = 5
          CALL VMAMSG (ISTAT,MAXFUN,ZERO)
          GO TO 999
        ENDIF
C
        NVMA = NVMA + 1
C
      ELSE
C
C...    CONTINUE TO CALCULATE NECESSARY INFORMATION FOR THE LINE SEARCH
C
        DIFFLS = FLS - FLSA
        DFLSA = DFLSAV
C
        IF (DIFFLS .LE. CHORD*DFLSA) THEN
          CALL GRADLF(G,M,N,IPRINT)
          BFLAG = .TRUE.
        ELSE
C
          IF (NVMA .LT. NVMLIN + NVMITR) THEN
            STEPL = ABIG(STPMAX,HALF*DFLSA/(DFLSA-DIFFLS))
C
C...        CALCULATE THE NEW CONTROL VARIABLE GUESSES AND SAVE THE OLD
C
            DO 40 I=1,N
              DELTA(I) = STEPL*DELTA(I)
              XA(I) = OXA(I)
              X(I) = XA(I) + DELTA(I)
   40       CONTINUE
C
            DFLSA = STEPL*DFLSA
            DFLSAV = DFLSA
C
            IF (NVMA .LT. MAXFUN) THEN
              NVMA = NVMA + 1
              GO TO 999
            ENDIF
          ENDIF
C
          DO 50 I=1,N
            X(I) = XA(I)
   50     CONTINUE
C
          NVMA = NVMA + 1
        ENDIF
      ENDIF
C
  999 CONTINUE
C
      IF (IPRINT .EQ. 2) THEN
        IBUG = 3
        CALL BUGAID (IBUG,IDUMMY,FLS)
        IBUG = 1
        CALL BUGAID (IBUG,IDUMMY,DFLSA)
        IF (.NOT. BFLAG) THEN
          IBUG = 4
          CALL BUGAID (IBUG,IDUMMY,STEPL)
          IBUG = 5
          CALL BUGAID (IBUG,N,XDUMMY)
        ENDIF
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE MATIN (C,NX,N,D,M,KEY,DETRMX )
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         DERIVE THE INVERSE RANK AND DETERMINANT OF AN INPUT
CD1         MATRIX (C). IF DESIRED THE INVERTED MATRIX CAN BE MULTIPLIED
CD1         BY A VECTOR (D) BY TURNING ON THE SWITCH M=1. THIS VERSION
CD1         CONVERTS SINGLE PRECISION INPUTS TO DOUBLE PRECISION ACCUR-
CD1         ACY INTERNALLY. WHEN MATRIX INVERSION IS COMPLETED THE
CD1         OUTPUTS ARE CONVERTED BACK TO SINGLE PRECISION.
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      C       THE MATRIX TO BE INVERTED
CD2      D       VECTOR TO MULTIPLY TIMES THE INVERTED MATRIX
CD2      DETERM  DETERMINANT OF THE MATRIX
CD2      KEY     RANK OF THE INVERTED MATRIX
CD2      M       SWITCH FOR THE VECTOR MULTIPLICATION
CD2      N       DIMENSION OF THE MATRIX A AND VECTOR B
CD2      NX      TOTAL SIZE OF THE MATRIX COMING IN
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      NONE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      NONE
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      TBD
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      NONE
C
C***********************************************************************
C
      DOUBLE PRECISION A
      DOUBLE PRECISION AMAX
      DOUBLE PRECISION B
      DOUBLE PRECISION DETERM
      DOUBLE PRECISION PIVOT
      DOUBLE PRECISION SWAP
      DOUBLE PRECISION T
      DOUBLE PRECISION ZERO
C
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXA   = MAXX * 2 )
C
      DIMENSION A      (MAXA  ,MAXA  )
      DIMENSION B      (     1)
      DIMENSION C      (NX    ,NX    )
      DIMENSION D      (     1)
      DIMENSION IPIVOT (MAXA  )
      DIMENSION INDEX  (MAXA  ,     2)
      DIMENSION PIVOT  (MAXA  )
C
C***********************************************************************
C
C...CONVERT THE C MATRIX AND B VECTOR TO DOUBLE PRECISON (A AND B)
C
      B(1) =  (D(1))
      DETERM =  (DETRMX)
C
      DO 15 I = 1,N
        DO 10 J=1,N
          A(I,J) =  (C(I,J))
   10   CONTINUE
   15 CONTINUE
C
C...INITIALIZATION
C
      ZERO = 1.D-16
      KEY = N
      DETERM = 1.D+0
      DO 20 J=1,N
        IPIVOT(J) = 0
   20 CONTINUE
C
      DO 550 I=1,N
C
C...    SEARCH FOR PIVOT ELEMENT
C
        AMAX = 0.D+0
C
        DO 105 J=1,N
          IF (IPIVOT(J) - 1 .NE. 0) THEN
            DO 100 K=1,N
              IF (IPIVOT(K) - 1 .LT. 0) THEN
                IF (DABS(AMAX) - DABS(A(J,K)) .LT. 0.) THEN
                  IROW = J
                  ICOLUM = K
                  AMAX = A(J,K)
                ENDIF
              ELSEIF (IPIVOT(J) - 1 .GT. 0) THEN
                GO TO 740
              ENDIF
  100       CONTINUE
          ENDIF
  105   CONTINUE
C
        IF (DABS(AMAX) .LE. ZERO ) THEN
          DETERM = 0.D0
          KEY = I-1
          GO TO 740
        ENDIF
        IPIVOT(ICOLUM) = IPIVOT(ICOLUM) + 1
C
C...    INTERCHANGE ROWSTO PUT PIVOT ELEMENT ON THE DIAGONAL
C
        IF (IROW - ICOLUM .NE. 0) THEN
          DETERM = -DETERM
          DO 200 L=1,N
            SWAP = A(IROW,L)
            A(IROW,L) = A(ICOLUM,L)
            A(ICOLUM,L) = SWAP
  200     CONTINUE
C
          IF (M .NE. 0) THEN
            SWAP = B(IROW)
            B(IROW) = B(ICOLUM)
            B(ICOLUM) = SWAP
          ENDIF
        ENDIF
C
        INDEX(I,1) = IROW
        INDEX(I,2) = ICOLUM
        PIVOT(I) = A(ICOLUM,ICOLUM)
C
C...    DIVIDE PIVOT ROW BY PIVOT ELEMENT
C
        A(ICOLUM,ICOLUM) = 1.D+0
C
        DO 350 L=1,N
          A(ICOLUM,L) = A(ICOLUM,L)/PIVOT(I)
  350   CONTINUE
C
        IF (M .NE. 0) THEN
          B(ICOLUM) = B(ICOLUM)/PIVOT(I)
        ENDIF
C
C...    REDUCE NON-PIVOT ROWS
C
        DO 540 L1=1,N
          IF (L1 - ICOLUM .NE. 0) THEN
            T = A(L1,ICOLUM)
            A(L1,ICOLUM) = 0.D0
            DO 450 L=1,N
              A(L1,L) = A(L1,L) - A(ICOLUM,L)*T
  450       CONTINUE
C
            IF (M .NE. 0) THEN
              B(L1) = B(L1) - B(ICOLUM)*T
            ENDIF
          ENDIF
  540   CONTINUE
  550 CONTINUE
C
C...INTERCHANGE COLUMNS
C
      DO 710 I=1,N
        L= N+ 1 - I
        IF (INDEX(L,1) - INDEX(L,2) .NE. 0) THEN
          JROW = INDEX (L,1)
          JCOLUM = INDEX(L,2)
          DO 705 K=1,N
            SWAP = A(K,JROW)
            A(K,JROW) = A(K,JCOLUM)
            A(K,JCOLUM) = SWAP
  705     CONTINUE
        ENDIF
  710 CONTINUE
C
      DO 720 I=1,N
        J = N + 1 - I
        DETERM = DETERM*PIVOT(J)
  720 CONTINUE
C
C...CONVERT BACK TO SINGLE PRECISION ACCURACY
C
  740 CONTINUE
      D(1) = (B(1))
      DETRMX = (DETERM)
      DO 770 I=1,N
        DO 760 J=1,N
          C(I,J) = (A(I,J))
  760   CONTINUE
  770 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE MOVE
CD0
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         STEPS ALONG THE SEARCH DIRECTION AND DISTANCE CALCULATED
CD1         TO GET A NEW VERTEX. IT ALSO CALCULATES NECESSARY INFORMA-
CD1         TION FOR ADDING OR EXCHANGING A CONSTRAINT IN THE BASIS
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD2
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      LSIZE   TEMPORARY ARRAY SIZE OF A VECTOR BEING MULTIPLIED
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      AMDAMX  THE MAXIMUM LAGRANGE MULTIPLIER
CD4      CAC     CURVATURE OF THE CONSTRAINT VECTOR CT (= CT*A*C)
CD4      CN      MATRIX OF CONSTRAINT NORMALS
CD4      CNEAR   DISTANCE TO THE NEAREST CURRENTLY INACTIVE CONSTRAINT
CD4      CSTAR   INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD4              CONSTRAINTS IN THE BASIS. SEE EQN (11)
CD4      CSTEP   COMPONENT OF THE GRADIENT ALONG THE SEARCH DIRECTION
CD4      H       CAN BE CONSIDERED AS REDUCED INVERSE HESSIAN OPERATOR
CD4              APPROPRIATE TO THE MANIFOLD FORMED BY THE INTERSECTION
CD4              OF THE CONSTRAINTS. SEE EQN (12)
CD4      ICNEAR  INDEX OF THE NEAREST CURRENTLY INACTIVE CONSTRAINT
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN
CD4              THE BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD4      NVAR2   TWICE THE VALUE OF NVAR. (= 2*(N+1))
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CCGCAC  REQUIRED DIRECTION OF SEARCH
CD5      CHC     PRODUCT OF THE NORMAL OF THE CONSTRAINT AND HDOTC
CD5      CSTARC  PRODUCT OF THE CSTAR OPERATOR AND NORMAL OF THE
CD5              CONSTRAINTS
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD5      HDOTC   PRODUCT OF THE H OPERATOR AND THE NORMAL OF THE
CD5              CONSTRAINTS
CD5      IACTIV  STATUS OF EACH OF THE MTOTAL CONSTRAINTS
CD5              =-1 EQUALITY
CD5              = 0 ACTIVE
CD5              = 1 INACTIVE
CD5              = 2 VIOLATED
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      MAXACT  MAXIMUM SIZE OF ARRY DEALING WITH THE IACTIVE ARRAY
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROLS AND THE
CD6              ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXACT = MAXX*2 + MAXC )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
C
      DIMENSION CSTAR  (MAXX  ,MAXC  )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMACT/ IACTIV (MAXACT)
C
      EQUIVALENCE ( CSTAR      ,  CINVRS ( 1,1) )
      EQUIVALENCE ( AMDAMX     ,  COMVMA (   1) )
      EQUIVALENCE ( CAC        ,  COMVMA (   2) )
      EQUIVALENCE ( CCGCAC     ,  COMVMA (   4) )
      EQUIVALENCE ( CHC        ,  COMVMA (   5) )
      EQUIVALENCE ( CNEAR      ,  COMVMA (   6) )
      EQUIVALENCE ( ICNEAR     ,  KOMVMA (   1) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
      EQUIVALENCE ( NVAR2      ,  KOMVMA (   7) )
C
C***********************************************************************
C
      DO 10 I=1,NVAR
        DELTA(I) = DELTA(I) + CNEAR*CSTEP(I)
   10 CONTINUE
C
      CCGCAC = CNEAR*AMDAMX/CAC
C
      IF (ICNEAR .LE. NVAR) THEN
        IF (NVAR .NE. NC) THEN
          DO 20 I=1,NVAR
            HDOTC(I) = H(I,ICNEAR)
   20     CONTINUE
        ENDIF
C
        DO 30 I=1,NC
          CSTARC(I) = CSTAR(I,ICNEAR)
   30   CONTINUE
        CHC = HDOTC(ICNEAR)
C
      ELSEIF (ICNEAR .LE. NVAR2) THEN
        J = ICNEAR - NVAR
        IF (NVAR .NE. NC) THEN
          DO 40 I=1,NVAR
            HDOTC(I) = -H(I,J)
   40     CONTINUE
        ENDIF
C
        DO 50 I=1,NC
          CSTARC(I) = -CSTAR(I,J)
   50   CONTINUE
        CHC = -HDOTC(J)
C
      ELSE
        J = ICNEAR - NVAR2
        LSIZE = 1 + MAXX*(NVAR - 1)
        IF (NVAR .NE. NC) THEN
          DO 60 I=1,NVAR
            HDOTC(I) = DOTPRD (NVAR, H(I,1),LSIZE,MAXX,CN(1,J),NVAR,1)
   60     CONTINUE
          CHC = DOTPRD (NVAR,CN(1,J),NVAR,1,HDOTC(1),NVAR,1)
        ENDIF
C
        DO 70 I=1,NC
          CSTARC(I) = DOTPRD (NVAR,CSTAR(I,1),LSIZE,MAXX,CN(1,J),NVAR,1)
   70   CONTINUE
      ENDIF
C
      IACTIV(ICNEAR) = 0
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE NEWBAS (ICROW)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         A NEW BASIS IS OBTAINED BY EXCHANGING CONSTRAINTS
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      ICROW   INDEX OF THE ROW OF CINVRS WHICH IS TO BE REMOVED FROM
CD2              THE BASIS
CD2
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      LSIZE   TEMPORARY ARRAY SIZE FOR A VECTOR BEING MULTIPLIED
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      CINVRS  INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD4              CONSTRAINTS IN THE BASIS. SEE EQN (2)
CD4      CN      MATRIX OF CONSTRAINT NORMALS
CD4      MTOTAL  TOTAL NUMBER OF CONSTRAINTS (= 2*(N+1)). INCLUDE THE
CD4              INPUT CONSTRAINTS AND THE UPPER AND LOWER BOUNDARIES
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD4      NVAR2   TWICE THE VALUE OF NVAR. (= 2*(N+1))
CD4      RESID   CONSTRAINT RESIDUAL
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD5      IACTIV  STATUS OF EACH OF THE MTOTAL CONSTRAINTS IN THE BASIS
CD5              =-1 EQUALITY
CD5              = 0 ACTIVE
CD5              = 1 INACTIVE
CD5              = 2 VIOLATED
CD5      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      ALPHA   NEAREST NONVIOLATED NONBASIC CONSTRAINT DISTANCE
CD6      BETA    FURTHEST VIOLATED CONSTRAINT DISTANCE AND THEN THE
CD6              MINIMUM BETWEEN ALPHA AND BETA
CD6      CROW    CONSTRAINT VECTOR OF CINVRS WHICH IS BEING REMOVED
CD6              FROM THE BASIS
CD6      DISTNC  DISTANCE MOVED AWAY FROM THE VERTEX IN THE DIRECTION
CD6              OF UI
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      IAINDX  INDEX OF THE CONSTRAINT CORRESPONDING TO ALPHA
CD6      IBINDX  INDEX OF THE CONSTRAINT CORRESPONDING TO BETA
CD6      MAXACT  MAXIMUM SIZE OF ARRAY DEALING WITH IACTIV ARRAY
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6              PLUS THE ARTIFICIAL VARIABLE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      RECIPR  RECIPRICAL OF THE RESIDUAL CONSTRAINT TO BE REMOVED
CD6              FROM THE BASIS
CD6      RRECIP  PRODUCT OF THE ITH RESIDUAL IN THE BASIS AND RECIPR
CD6      UI      ITH ROW OF CINVRS. DIRECTION OF SEARCH FROM THE
CD6              CURRENT VERTEX
CD6      UICJ    PRODUCT OF THE ITH ROW OF CINVRS AND THE NORMAL VECTORS
CD6              OF CONSTRAINTS IN THE BASIS
CD6      XINF    NUMBER FOR INFINITY
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "THE CALCULATION OF FEASIBLE
CD7      POINTS FOR LINEARLY CONSTRAINED OPTIMIZATION PROBLEMS" AERE-
CD7      R.6354 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXACT = MAXX*2 + MAXC )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( XINF   = 1.E38 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION CROW   (MAXX  )
      DIMENSION UI     (MAXX  )
C
      COMMON   /CVMACT/ RESID  (MAXACT)
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMACT/ IACTIV (MAXACT)
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( MTOTAL     ,  KOMVMA (   4) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
      EQUIVALENCE ( NVAR2      ,  KOMVMA (   7) )
C
C***********************************************************************
C
C...MOVE AWAY FROM THE VERTEX IN THE DIRECTION OF UI
C
      DO 10 I=1,NVAR
        UI(I) = CINVRS(ICROW,I)
   10 CONTINUE
C
      ALPHA = XINF
      BETA = ZERO
C
      DO 20 I=1,MTOTAL
        IF (IACTIV(I) .GE. 1) THEN
          IF (I .LE. NVAR) THEN
            UICJ = -UI(I)
          ELSEIF (I .LE. NVAR2) THEN
            J = I - NVAR
            UICJ = UI(J)
          ELSE
            J = I - NVAR2
            UICJ = -DOTPRD (NVAR,UI(1),NVAR,1,CN(1,J),NVAR,1)
          ENDIF
C
          IF (IACTIV(I) .EQ. 1) THEN
C
C...        FIND THE NEAREST NONVIOLATED NONBASIS CONSTRAINT
C
            IF (UICJ .GT. 0.) THEN
              DISTNC = RESID(I)/UICJ
              IF (DISTNC .LT. ALPHA) THEN
                ALPHA = DISTNC
                IAINDX = I
              ENDIF
            ENDIF
          ELSE
C
C...        FIND THE FURTHEST VIOLATED CONSTRAINT
C
            IACTIV(I) = 1
C
            IF (UICJ .LT. 0.) THEN
              DISTNC = RESID(I)/UICJ
              IF (DISTNC .GT. BETA) THEN
                BETA = DISTNC
                IBINDX = I
              ENDIF
            ENDIF
          ENDIF
        ENDIF
   20 CONTINUE
C
      IF (ALPHA .LE. BETA) THEN
C
C...    NEW VERTEX AT X + ALPHA*UI
C
        IBINDX = IAINDX
        BETA = ALPHA
      ENDIF
C
C...EXCHANGE THE CONSTRAINTS WITHIN THE BASIS
C
      IACTIV(IBASIS(ICROW)) = 1
      IACTIV(IBINDX) = 0
      IBASIS(ICROW) = IBINDX
C
C...UTILIZE THE SIMPLEX METHOD TO EXCHANGE WITHIN CINVRS
C
      IF (IBINDX .LE. NVAR) THEN
        DO 30 I=1,NVAR
          CROW(I) = CINVRS(I,IBINDX)
   30   CONTINUE
      ELSEIF (IBINDX .LE. NVAR2) THEN
        IBINDX = IBINDX - NVAR
        DO 40 I=1,NVAR
          CROW(I) = -CINVRS(I,IBINDX)
   40   CONTINUE
      ELSE
        IBINDX = IBINDX - NVAR2
        LSIZE = 1 + MAXX*(NVAR - 1)
        DO 50 I=1,NVAR
          CROW(I) = DOTPRD (NVAR,CINVRS(I,1),LSIZE,MAXX,CN(1,IBINDX),
     1                      NVAR,1)
   50   CONTINUE
      ENDIF
C
      RECIPR = FP1/CROW(ICROW)
C
      DO 80 I=1,NVAR
        DELTA(I) = DELTA(I) + BETA*UI(I)
        IF (I .EQ. ICROW) THEN
          DO 60 J=1,NVAR
            CINVRS(I,J) = CINVRS(I,J)*RECIPR
   60     CONTINUE
        ELSE
          RRECIP = RECIPR*CROW(I)
          DO 70 J=1,NVAR
            CINVRS(I,J) = CINVRS(I,J) - RRECIP*UI(J)
   70     CONTINUE
        ENDIF
   80 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE NEWDIR
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         CALCULATES THE NEW MAXIMUM LAGRANGE MULTIPLIER, GRADIENT
CD1         DIRECTION, AND SEARCH DIRECTION
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      LSIZE   TEMPORARY ARRAY SIZE FOR A VECTOR BEING MULTIPLIED
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      A       VARIABLE METRIC MATRIX. CORRESPONDS TO THE HESSIAN
CD4              MATRIX IN THE QUADRATIC FUNCTION. SEE EQN (1A)
CD4      CSTEP   COMPONENT OF THE GRADIENT IN THE SEARCH DIRECTION
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      GM      NEGATIVE OF THE GRADIENT VECTOR ( =-G). CORRESPONDS
CD4              TO B IN EQNS  (4),(5),(6), AND (8)
CD4      H       CAN BE CONSIDERED AS A REDUCED INVERSE HESSIAN OPERATOR
CD4              APPROPRIATE TO THE MANIFORLD FORMED BY THE INTERSECTION
CD4              OF THE CONSTRAINTS. SEE EQN (12)
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      AMDAMX  THE MAXIMUM LAGRANGE MULTIPLIER
CD5      GHAT    CURRENT GRADIENT VECTOR OF THE OBJECTIVE FUNCTION.
CD5              SEE EQN (9)
CD5      PASSIV  INDICATES A CONSTRAINT HAS BEEN REMOVED TO GET A
CD5              SEARCH DIRECTION
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS AND
CD6              THE ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      TOL     TOLERANCE FOR EQUALITY OF ZERO
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      LOGICAL PASSIV
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( TOL    = 1.E-38 )
      PARAMETER ( ZERO   =  0.0 )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( AMDAMX     ,  COMVMA (   1) )
      EQUIVALENCE ( PASSIV     ,  COMVMA (   8) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
      LSIZE = 1 + MAXX*(NVAR - 1)
      DO 10 I=1,NVAR
        GHAT(I) = DOTPRD(NVAR,A(I,1),LSIZE,MAXX,DELTA(1),NVAR,1) - GM(I)
   10 CONTINUE
C
      AMDAMX = ZERO
C
      DO 20 I=1,NVAR
        CSTEP(I) = -DOTPRD (NVAR,H(I,1),LSIZE,MAXX,GHAT(I),NVAR,1)
        IF (ABS(CSTEP(I)) .GT. TOL) THEN
          AMDAMX = FP1
        ENDIF
   20 CONTINUE
C
      PASSIV = .FALSE.
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE NODESG
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         NO DESIGNATED CONSTRAINTS, VERTEX CHOSEN FROM THE LOWER
CD1         AND UPPER BOUNDS. (WHICHEVER IS CLOSER)
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD2
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      BDL     LOWER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1A)
CD4      BDU     UPPER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1A)
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD5      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CINVRS  INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD5              CONSTRAINTS IN THE BASIS. SEE EQN (2)
CD5      IACTIV  STATUS OF EACH OF THE MTOTAL CONSTRAINTS
CD5              =-1 EQUALITY
CD5              = 0 ACTIVE
CD5              = 1 INACTIVE
CD5              = 2 VIOLATED
CD5      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD5      NC      NUMBER OF CONSTRAITS IN THE BASIS. INITIALLY
CD5              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN
CD5              THE BASIS.
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXACT  MAXIMUM SIZE OF ARRAY DEALING WITH THE IACTIV ARRAY
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROLS AND THE
CD6              ARTIFICIAL VARIABLE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "THE CALCULATION OF FEASIBLE
CD7      POINTS FOR LINEARLY CONSTRAINED OPTIMIZATION PROBLEMS" AERE-
CD7      R.6354 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXACT = MAXX*2 + MAXC )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMACT/ IACTIV (MAXACT)
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
      DO 20 I=1,NVAR
        DO 10 J=1,NVAR
          CINVRS(I,J) = ZERO
   10   CONTINUE
C
        IF (DELTA(I) - BDL(I) .LE. BDU(I) - DELTA(I)) THEN
          IBASIS(I) = I
          CINVRS(I,I) = FP1
        ELSE
          IBASIS(I) = NVAR + I
          CINVRS(I,I) = -FP1
        ENDIF
C
        J = IBASIS(I)
        IACTIV(J) = 0
   20 CONTINUE
C
      NC = NVAR
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE QUAD (ISTAT,MEQ,IPRINT)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         MINIMIZES A QUADRATIC FUNCTION OF N VARIABLES SUBJECT TO M
CD1         LINEAR EQUALITY AND INEQUALITY CONSTRAINTS.
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      HCMAT   GENERAL MATRIX CONTAINING EQUATIONS GOVERNING EQUALITY
CD2              PROBLEM EQN (7). PARTITIONS OF THIS MATRIX ARE THE H
CD2              AND CSTAR OPERATORS.
CD2      IPRINT  FLAG TO INDICATE AMOUNT OF PRINTOUT
CD2              = 0 NO PRINTOUT
CD2              = 1 VALUES OF X, C, AND F ARE PRINTED
CD2              = 2 DEBUGGING PRINTOUT
CD2      ISTAT   INDEX TO INDICATE THE STATUS OF THE OPTIMIZATION
CD2              ALGORITHM
CD2              =-2 INITIAL PASS THROUGH VMACO WITH THE B MATRIX
CD2                  ALREADY INITIALIZED
CD2              =-1 INITIAL PASS THROUGH VMACO
CD2              = 0 CONTINUE ALGORITH CALCULATIONS
CD2              = 1 CONVERGE WITHIN REQUIRED ACCURACY
CD2              > 1 ERROR ENCOUNTERED
CD2      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD2      RESIDU  CONSTRAIN RESIDUAL FOR THE QUADRATIC EQUALITY PROBLEM
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      AMDAMX  THE MAXIMUM LAGRANGE MULTIPLIER
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN
CD4              THE BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CSTAR   INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD5              CONSTRAINTS IN THE BASIS. SEE EQN (11)
CD5      H       CAN BE CONSIDERED AS A REDUCED INVERSE HESSIAN OPERATOR
CD5              APPROPRIATE TO THE MANIFOLD FORMED BY THE INTERSECTION
CD5              OF THE CONSTRAINTS. SEE EQN (12)
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      MAXHC   MAXIMUM SIZE OF SQUARE ARRAY DEALING WITH HCMAT
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6              AND THE ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      NSERCH  COUNTER ON THE ITERATION NUMBER CURENTLY ON IN THE
CD6              SEARCH FOR THE MINIMUM OF THE QUADRATIC FUNCTION
CD6      NSCHMX  MAXIMUM VALUE ALLOWED FOR NSERCH
CD6      RETEST  FLAG TO INDICATE WHETHER THE FINAL SOLUTION HAS BEEN
CD6              VERIFIED
CD6      TOL     TOLERANCE FOR EQUALITY OF ZERO
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      LOGICAL RETEST
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXHC  = MAXX*2 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( NSCHMX =  100 )
      PARAMETER ( TOL    = 1.E-38 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION CSTAR  (MAXX  ,MAXC  )
      DIMENSION HCMAT  (MAXHC ,MAXHC )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( CSTAR      ,  CINVRS ( 1,1) )
      EQUIVALENCE ( AMDAMX     ,  COMVMA (   1) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
      RETEST = .FALSE.
C
C...CALL THE FEASIBLE VERTEX ROUTINE
C
   10 CONTINUE
      CALL VERTEX (MEQ)
C
      IF (ISTAT .NE. 0) THEN
        NC = 0
        GO TO 999
      ENDIF
C
C...INITIALIZE THE OPERATOR H=0.
C
      NC = NVAR
      DO 30 I=1,NVAR
        DO 20 J=1,NVAR
          H(I,J) = ZERO
   20   CONTINUE
   30 CONTINUE
C
      NSERCH = 0
   40 CONTINUE
      NSERCH = NSERCH + 1
C
C...PROTECT AGAINST AN INFINITE LOOP OCCURRING
C
      IF (NSERCH .GE. NSCHMX) THEN
        ISTAT = 14
        CALL VMAMSG (ISTAT,NSCHMX,XDUMMY)
        IBUG = 5
        CALL BUGAID (IBUG,NVAR,XDUMMY)
        IBUG = 7
        CALL BUGAID (IBUG,NVAR,XDUMMY)
        IBUG = 6
        CALL BUGAID (IBUG,IDUMMY,AMDAMX)
        GO TO 999
      ENDIF
C
C...INITIALIZE THE LAGRANGIAN FUNCTION
C
      CALL LGRANG
C
      IF (IPRINT .EQ. 2) THEN
        IBUG = 5
        CALL BUGAID (IBUG,NVAR,XDUMMY)
        IBUG = 7
        CALL BUGAID (IBUG,NVAR,XDUMMY)
        IBUG = 6
        CALL BUGAID (IBUG,IDUMMY,AMDAMX)
      ENDIF
C
      IF (AMDAMX .GT. ZERO) THEN
C
C...    DOWNHILL PROGRESS CAN BE MADE BY DROPPING A CONSTRAINT FROM
C...    THE BASIS
C
        CALL SERCH
        GO TO 40
      ENDIF
C
      IF (RETEST) THEN
C
C...    A VERIFIED SOLUTION HAS BEEN FOUND
C
        IF (ABS(AMDAMX) .LT. TOL) THEN
C
C...      SOLUTION MAY BE A DEGENERATE LOCAL MINIMUM
C
          CALL VMAMSG (10,0,ZERO)
        ENDIF
      ELSE
C
C...    RETEST THE LAGRANGE FUNCTION BECAUSE THE SUCCESSIVE UPDATING
C...    OF THE OPERATORS (CSTAR AND H) MIGHT CAUSE AN ACCUMULATION ERROR
C
        RETEST = .TRUE.
C
        CALL RESET (HCMAT,RESIDU)
C
        IF (RESIDU .LT. ZERO) THEN
          GO TO 10
        ELSE
          DO 60 I=1,NVAR
            IX = I+NVAR
            DO 50 J=1,NVAR
              H(I,J)=HCMAT(I,J)
              CSTAR(I,J)=HCMAT(IX,J)
   50       CONTINUE
   60     CONTINUE
          GO TO 40
        ENDIF
      ENDIF
C
  999 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE REMOVE
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         CONSTRAINT IS REMOVED FROM THE BASIS AND THE CORRESPONDING
CD1         H AND CSTAR OPERATORS ARE UPDATED
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD2
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      LSIZE   TEMPORARY ARRAY SIZE OF A VECTOR BEING MULTIPLIED
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      AC      PRODUCT OF THE A MATRIX CT VECTOR
CD4      CAC     CURVATURE OF THE CONSTRAINT VECTOR CT (= CT*A*C)
CD4      CT      CONSTRAINT VECTOR, CORRESPONDING TO THE LARGEST
CD4              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD4              BASIS
CD4      LAMROW  INDEX OF THE CONSTRAINT, CORRESPONDING TO THE LARGEST
CD4              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD4              BASIS
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN
CD4              THE BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CKAC    PRODUCT OF THE KTH ROW OF CSTAR MATRIX AND AC MATRIX
CD5      CSTAR   INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD5              CONSTRAINTS IN THE BASIS. SEE EQN (11)
CD5      H       CAN BE CONSIDERED AS A REDUCED INVERSE HESSIAN OPERATOR
CD5              APPROPRIATE TO THE MANIFOLD FORMED BY THE INTERSECTION
CD5              OF THE CONSTRAINTS. SEE EQN (12)
CD5      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      CKACDC  VALUE OF CKAC DIVIDED BY CAC
CD6      CTDCAC  VALUE OF CT DIVIDED BY CAC
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS AND
CD6              ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
C
      DIMENSION CSTAR  (MAXX  ,MAXC  )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( CSTAR      ,  CINVRS ( 1,1) )
      EQUIVALENCE ( CAC        ,  COMVMA (   2) )
      EQUIVALENCE ( LAMROW     ,  KOMVMA (   3) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
C...FIRST THE H OPERATOR IS UPDATED
C
      DO 20 I=1,NVAR
        CTDCAC = CT(I)/CAC
        DO 10 J=1,I
          H(I,J) = H(I,J) + CTDCAC*CT(J)
          H(J,I) = H(I,J)
   10   CONTINUE
   20 CONTINUE
C
      IF (NC .GT. 1) THEN
        IF (LAMROW .NE. NC) THEN
C
C...      REPLACE THE DELEGATED CONSTRAINT BY MOVING THE LAST CONSTRAINT
C...      IN THE BASIS INTO ITS PLACE
C
          DO 30 I=1,NVAR
            CSTAR(LAMROW,I) = CSTAR(NC,I)
   30     CONTINUE
          IBASIS(LAMROW) = IBASIS(NC)
        ENDIF
C
C...    UPDATE THE CSTAR OPERATOR MATRIX
C
        NC = NC - 1
C
        LSIZE = 1 + MAXX*(NVAR - 1)
        DO 40 I=1,NC
          CKAC(I) = DOTPRD (NVAR,CSTAR(I,1),LSIZE,MAXX,AC(1),NVAR,1)
   40   CONTINUE
C
        DO 60 I=1,NC
          CKACDC = CKAC(I)/CAC
          DO 50 J=1,NVAR
            CSTAR(I,J) =CSTAR(I,J) - CKACDC*CT(J)
   50     CONTINUE
   60   CONTINUE
C
      ELSE
        NC = 0
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE RESET (HCMAT,RESIDU)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         SETS UP THE MATRIX AND RIGHT HAND SIDE OF EQUATIONS
CD1         GOVERNING THE QUADRATIC (AND LAGRANGIAN FUNCTION) EQUALITY
CD1         PROBLEM
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      HCMAT   GENERAL MATRIX CONTAINING EQUATIONS GOVERNING EQUALITY
CD3              PROBLEM. SEE EQN (7). PARTITIONS OF THES MATRIX ARE
CD3              H AND CSTAR OPERATORS
CD3      KEY     RANK OF THE INVERTED MATRIX
CD3      RESIDU  CONSTRAINT RESIDUAL FOR THE QUADRATIC EQUALITY PROBLEM
CD3      VECDUM  DUMMY VALUE FOR THE ROUTINE MATIN
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      A       VARIABLE METRIC MATRIX. CORRESPONDS TO THE HESSIAN
CD4              MATRIX IN THE QUADRATIC FUNCTION. SEE EQN (1A)
CD4      BDL     LOWER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1B)
CD4      BDU     UPPER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1B)
CD4      CN      MATRIX OF CONSTRAINT NORMALS
CD4      D       RIGHT HAND SIDE OF CONSTRAINT EQUATION. SEE EQN (2)
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      GM      NEGATIVE OF THE GRADIENT VECTOR  ( =-G). CORRESPONDS
CD4              TO B IN EQNS (4),(5),(6), AND (8)
CD4      IACTIV  STATUS OF EACH OF THE MTOTAL CONSTRAINTS
CD4              =-1 EQUALITY
CD4              = 0 ACTIVE
CD4              = 1 INACTIVE
CD4              > 1 ERROR ENCOUNTERED
CD4      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD4      MTOTAL  TOTAL NUMBER OF CONSTRAINTS (= M+2*NVAR). INCLUDES THE
CD4              INPUT CONSTRAINTS AND THE UPPER AND LOWER BOUNDARIES
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN
CD4              THE BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD4      NVAR2   TWICE THE VALUE OF NVAR. (= 2*(N+1))
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      NONE
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXACT  MAXIMUM SIZE OF ARRAY DEALING WITH THE IACTIV ARRAY
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXHC   MAXIMUM SIZE OF ARRAYS DEALING WITH THE HCMAT MATRIX
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS PLUS
CD6              THE ARTIFICIAL VARIABLE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      NNC     SUM OF NVAR + NC
CD6      RHSEP   RIGHT HAND SIDE OF THE EUATION GOVERNING THE EQUALITY
CD6              PROBLEM
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXACT = MAXX*2 + MAXC )
      PARAMETER ( MAXHC  = MAXX*2 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION RHSEP  (MAXHC )
      DIMENSION HCMAT  (MAXHC ,MAXHC )
C
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMACT/ IACTIV (MAXACT)
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( MTOTAL     ,  KOMVMA (   4) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR        , KOMVMA (   6) )
      EQUIVALENCE ( NVAR2       , KOMVMA (   7) )
C
C***********************************************************************
C
      DO 20 I=1,NVAR
        RHSEP(I) = GM(I)
        DO 10 J=1,NVAR
          HCMAT(I,J) = A(I,J)
   10   CONTINUE
   20 CONTINUE
C
      IF (NC .NE. 0) THEN
        DO 60 I=1,NC
          IF (IBASIS(I) .LE. NVAR2) THEN
            DO 30 J=1,NVAR
              HCMAT(J,NVAR+I) = ZERO
              HCMAT(NVAR+I,J) = ZERO
   30       CONTINUE
C
            IF (IBASIS(I) .LE. NVAR) THEN
              HCMAT(NVAR+I,IBASIS(I)) = FP1
              HCMAT(IBASIS(I),NVAR+I) = FP1
              RHSEP(NVAR+I) = BDL(IBASIS(I))
            ELSE
              J = IBASIS(I) - NVAR
              HCMAT(NVAR+I,J) = -FP1
              HCMAT(J,NVAR+I) = -FP1
              RHSEP(NVAR+I) = -BDU(J)
            ENDIF
          ELSE
            J = IBASIS(I) - NVAR2
            DO 40 JX = 1,NVAR
              HCMAT(NVAR+I,JX) = CN(JX,J)
              HCMAT(JX,NVAR+I) = CN(JX,J)
   40       CONTINUE
            RHSEP(NVAR+I) = D(J)
          ENDIF
C
          DO 50 J=1,NC
            IX = NVAR + I
            JX = NVAR + J
            HCMAT(IX,JX) = ZERO
   50     CONTINUE
   60   CONTINUE
      ENDIF
C
      NNC = NVAR + NC
C
C...INSERT THE GENERAL MATRIX
C
      CALL MATIN (HCMAT,MAXHC,NNC,VECDUM,0,KEY,DETERM)
      IF (KEY .LT. NNC) THEN
        CALL VMAMSG (15,2,DUMMY)
      ENDIF
C
      DO 70 I=1,NVAR
        DELTA(I) = DOTPRD (NNC,HCMAT(1,I),NNC,1,RHSEP(1),NNC,1)
   70 CONTINUE
C
C...CHECK THE FEASIBILITY OF THE NEW POINT
C
      DO 80 I=1,MTOTAL
        IF (IACTIV(I) .GE. 1 ) THEN
          IF (I .LE. NVAR) THEN
            RESIDU = DELTA(I) - BDL(I)
          ELSEIF (I .LE. NVAR2) THEN
            J = I-NVAR
            RESIDU = BDU(J) - DELTA(J)
          ELSE
            J = I-NVAR2
            RESIDU = DOTPRD (NVAR,CN(1,J),NVAR,1,DELTA(1),NVAR,1) - D(J)
          ENDIF
        ENDIF
   80 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE SERCH
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         SEARCHES FOR AN OPTIMAL SOLUTION BY DROPPING A CONSTRAINT
CD1         FROM THE BASIS
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD2
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      DEPEND  FLAG TO INDICATE THE DEPENDENCY OF THE NEW CONSTRAINT
CD3      POSTIV  CURVATURE OF THE CONSTRAINT VECTOR BEING REMOVED FROM
CD3              THE BASIS
CD3              = .TRUE. POSITIVE CURVATURE (OPTIMAL)
CD3              = .FALSE. NEGATIVE CURVATURE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      AMDAMX  THE MAXIMUM LAGRANGE MULTIPLIER
CD4      CNEAR   DISTANCE TO THE NEAREST CURRENTLY INACTIVE CONSTRAINT
CD4      CSTAR   INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD4              CONSTRAINTS IN THE BASIS. SEE EQN (11)
CD4      CSTEP   COMPONENT OF THE GRADIENT ALONG THE SEARCH DIRECTION
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN THE
CD4              BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD4      PASSIV  INDICATES A CONSTRAINT HAS BEEN REMOVED TO GET A
CD4              SEARCH DIRECTION
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CT      CONSTRAINT VECTOR, CORRESPONDING TO THE LARGEST
CD5              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD5              BASIS
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD5      LAMROW  INDEX OF THE CONSTRAINT, CORRESPONDING TO THE LARGEST
CD5              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD5              BASIS
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS AND
CD6              ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      TOL     TOLERANCE EQUALITY FOR ZERO
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      LOGICAL DEPEND
      LOGICAL PASSIV
      LOGICAL POSTIV
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( TOL    = 1.E-38 )
C
      DIMENSION CSTAR  (MAXX  ,MAXC  )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( CSTAR      ,  CINVRS ( 1,1) )
      EQUIVALENCE ( AMDAMX     ,  COMVMA (   1) )
      EQUIVALENCE ( CNEAR      ,  COMVMA (   6) )
      EQUIVALENCE ( PASSIV     ,  COMVMA (   8) )
      EQUIVALENCE ( LAMROW     ,  KOMVMA (   3) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
C...SET THE DIRECTION OF SEARCH AS CORRESPONDING TO THE ROW OF CSTAR
C...WITH THE LARGEST POSITIVE LAGRANGE MULTIPLIER
C
      DO 10 I=1,NVAR
        CT(I) = CSTAR(LAMROW,I)
   10 CONTINUE
C
   20 CONTINUE
      CALL STEPC (POSTIV)
C
   30 CONTINUE
      CALL SLNEAR
C
      IF (POSTIV .AND. CNEAR .GE. FP1) THEN
C
C...    A MINIMUM HAS BEEN FOUND
C
        DO 40 I=1,NVAR
          DELTA(I) = DELTA(I) + CSTEP(I)
   40   CONTINUE
        IF (PASSIV) THEN
          CALL REMOVE
        ENDIF
        GO TO 999
      ENDIF
C
C...MOVE TO A NEW VERTEX AND CALCULATE PERTINENT INFORMATION
C
      CALL MOVE
C
      IF (NC .EQ. NVAR) THEN
C
C...    APPLY THE SIMPLEX FORMULA TO EXCHANGE CONSTRAINTS
C
        CALL EXCHNG
        GO TO 999
      ENDIF
C
      DEPEND = .TRUE.
      IF (PASSIV) THEN
C
C...    CHECK THE DEPENDENCY OF THE NEW CONSTRAINT
C
        CALL CHKDEP (DEPEND)
      ENDIF
C
      IF (DEPEND) THEN
C
C...    APPLY FORMULA FOR EXCHANGING NEW CONSTRAINT WITH A PASSIVE ONE
C
        CALL CHANGP
      ELSE
C
C...    APPLY FORMULA FOR ADDING A CONSTRAINT TO THE BASIS
C
        CALL ADD
        IF (PASSIV) THEN
C
C...      REMOVAL OF A CONSTRAINT HAS BEEN DEFERRED. SET UP AS IF THE
C...      CONSTRAINT IS BEING REMOVED FROM THE AUGMENTED BASIS
C
          CALL DEFER
C
          IF (ABS(AMDAMX) .LT. TOL) THEN
            CALL REMOVE
            GO TO 999
          ELSE
            GO TO 20
          ENDIF
        ENDIF
      ENDIF
C
      IF (NC .NE. NVAR) THEN
C
C...    CALCULATE THE NEW SEARCH DIRECTION
C
        CALL NEWDIR
        IF (ABS(AMDAMX) .GT. TOL) THEN
          POSTIV = .TRUE.
          GO TO 30
        ENDIF
      ENDIF
C
  999 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE SLNEAR
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         A LINEAR SEARCH ALONG THE DIRECTION OF SEARCH IS CONDUCTED
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      BDL     LOWER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1B)
CD4      BDU     UPPER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1B)
CD4      CN      MATRIX OF CONSTRAINT NORMALS
CD4      CSTEP   COMPONENT OF THE GRADIENT ALONG THE SEARCH DIRECTION
CD4      D       RIGHT HAND SIDE OF CONSTRAINT EQUATION. SEE EQN (2)
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD4      LAMROW  INDEX OF THE CONSTRAINT, CORRESPONDING TO THE LARGEST
CD4              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD4              BASIS
CD4      MTOTAL  TOTAL NUMBER OF CONSTRAINTS (= M+2*NVAR). INCLUDES THE
CD4              INPUT CONSTRAINTS AND THE UPPER AND LOWER BOUNDARIES
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIBLE
CD4      NVAR2   TWICE THE VALUE OF NVAR. (= 2*(N+1))
CD4      PASSIV  INDICATES A CONSTRAINT HAS BEEN REMOVED TO GET A
CD4              SEARCH DIRECTION
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      CNEAR   DISTANCE TO THE NEAREST CURRENTLY INACTIVE CONSTRAINT
CD5      IACTIV  STATUS OF EACH OF THE MTOTAL CONSTRAINTS
CD5              =-1 EQUALITY
CD5              = 0 ACTIVE
CD5              = 1 INACTIVE
CD5              = 2 VIOLATED
CD5      ICNEAR  INDEX OF THE NEAREST CURRENTLY INACTIVE CONSTRAINT
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      MAXACT  MAXIMUM SIZE OF ARRAY DEALING WITH THE IACTIVE ARRAY
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROLS PLUS THE
CD6              ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      XINF    NUMBER FOR INFINITY
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      LOGICAL PASSIV
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXACT = MAXX*2 + MAXC )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( XINF   = 1.E38)
      PARAMETER ( ZERO   =  0.0 )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMACT/ IACTIV (MAXACT)
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( CNEAR      ,  COMVMA (   6) )
      EQUIVALENCE ( PASSIV     ,  COMVMA (   8) )
      EQUIVALENCE ( ICNEAR     ,  KOMVMA (   1) )
      EQUIVALENCE ( LAMROW     ,  KOMVMA (   3) )
      EQUIVALENCE ( MTOTAL     ,  KOMVMA (   4) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
      EQUIVALENCE ( NVAR2      ,  KOMVMA (   7) )
C
C***********************************************************************
C
      CNEAR = XINF
C
      DO 10 I=1,MTOTAL
        IF (IACTIV(I) .EQ. 1) THEN
C
C...      THE CONSTRAINT BEING CONSIDERED IS INACTIVE
C
          IF (I .LE. NVAR) THEN
            IF (CSTEP(I) .LT. ZERO) THEN
              DISTC = (BDL(I)-DELTA(I))/CSTEP(I)
              IF (DISTC .LT. CNEAR) THEN
                CNEAR = DISTC
                ICNEAR = I
              ENDIF
            ENDIF
          ELSEIF (I .LE. NVAR2) THEN
            J = I - NVAR
            IF (CSTEP(J) .GT. ZERO) THEN
              DISTC = (BDU(J) - DELTA(J))/CSTEP(J)
              IF (DISTC .LT. CNEAR) THEN
                CNEAR = DISTC
                ICNEAR = I
              ENDIF
            ENDIF
          ELSE
            J = I - NVAR2
            CCSTEP = DOTPRD (NVAR,CN(1,J),NVAR,1,CSTEP(1),NVAR,1)
            IF (CCSTEP .LT. ZERO) THEN
              DISTC = D(J) -DOTPRD (NVAR,CN(1,J),NVAR,1,DELTA(1),NVAR,1)
              DISTC = DISTC/CCSTEP
              IF (DISTC .LT. CNEAR) THEN
                CNEAR = DISTC
                ICNEAR = I
              ENDIF
            ENDIF
          ENDIF
        ENDIF
   10 CONTINUE
C
C...THE CURRENT CONSTRAINT (DIRECTION OF SEARCH) IS REMOVED FROM THE
C...BASIS
C
      IF (PASSIV) THEN
        IACTIV(IBASIS(LAMROW)) = 1
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE STEPC (POSTIV)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         CALCULATES THE STEP TO THE STATIONARY POINT IN THE NEW BASIS
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      NONE
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      LSIZE   TEMPORARY ARRAY SIZE OF A VECTOR BEING MULTIPLIED
CD3      POSTIV  CURVATURE OF THE CONSTRAINT VECTOR BEING REMOVED FROM
CD3              THE BASIS
CD3              = .TRUE. POSITIVE CURVATURE (OPTIMAL)
CD3              = .FALSE. NEGATIVE CURVATURE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      A       VARIABLE METRIC MATRIX. CORRESPONDS TO THE HESSIAN
CD4              MATRIX IN THE QUADRATIC FUNCTION. SEE EQN (1A)
CD4              EQUIVALENT TO THE B MATRIX IN VMACO
CD4      AMDAMX  THE MAXIMUM LAGRANG MULTIPLIER
CD4      CT      CONSTRAINT VECTOR, CORRESPONDING TO THE LARGEST
CD4              LAGRANGE MULTIPLIER, WHICH IS TO BE REMOVED FROM THE
CD4              BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      AC      PRODUCT OF THE A MATRIX AND THE CT VECTOR
CD5      CAC     CURVATURE OF THE CONSTRAINT VECTOR CT (= CT*A*C)
CD5      CSTEP   COMPONENT OF THE GRADIENT ALONG THE SEARCH DIRECTION
CD5      PASSIV  INDICATES A CONSTRAINT HAS BEEN REMOVED TO GET A
CD5              SEARCH DIRECTION
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROLS PLUS THE
CD6              ARTIFICIAL VARIABLE
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITH CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      LOGICAL PASSIV
      LOGICAL POSTIV
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX3/ AC     (MAXX  )
     1   ,              CKAC   (MAXX  )
     1   ,              CSTARC (MAXX  )
     1   ,              CSTEP  (MAXX  )
     1   ,              CT     (MAXX  )
     1   ,              HDOTC  (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( AMDAMX     ,  COMVMA (   1) )
      EQUIVALENCE ( CAC        ,  COMVMA (   2) )
      EQUIVALENCE ( PASSIV     ,  COMVMA (   8) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
C
C***********************************************************************
C
      LSIZE = 1 + MAXX*(NVAR - 1)
      DO 10 I=1,NVAR
        AC(I) = DOTPRD (NVAR,A(I,1),LSIZE,MAXX,CT(1),NVAR,1)
   10 CONTINUE
      CAC = DOTPRD (NVAR,CT(1),NVAR,1,AC(1),NVAR,1)
C
      IF (CAC .LE. ZERO) THEN
C
C...    CURVATURE OF THE CONSTRAINT VECTOR BEING REMOVE FROM THE BASIS
C...    IS NEGATIVE AND THEREFORE THE STEP TO A STATIONARY POINT IS
C...    NOT A MINIMUM
C
        POSTIV = .FALSE.
        DO 20 I=1,NVAR
          CSTEP(I) = CT(I)
   20   CONTINUE
      ELSE
C
C...    CURVATURE OF CONSTRAINT VECTOR BEING REMOVED FROM THE BASIS
C...    IS POSITIVE AND THE MINIMUM CAN BE FOUND AT CSTEP
C
        POSTIV = .TRUE.
        DO 30 I=1,NVAR
          CSTEP(I) = CT(I)*AMDAMX/CAC
   30   CONTINUE
      ENDIF
C
      PASSIV = .TRUE.
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE STINIT (ISTAT,M,MEQ,N,VLARGE)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         CURRENTLY ON THE INITIAL CALL TO VMACO (ISTAT=-1) AND
CD1         THEREFORE VARIABLES AND ELEMENTS ARE SET UP ACCORDINGLY
CD1         FOR THE QUADRATIC PROGRAMMING FUNCTION
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      ISTAT   INDEX TO INDICATE THE STATUS OF THE OPTIMIZATION
CD2              ALGORITHM
CD2              =-2 INITIAL PASS THROUGH VMACO WITH THE B MATRIX
CD2                  ALREADY INITIALIZED
CD2              =-1 INITIAL PASS THROUGH VMACO
CD2              = 0 CONTINUE ALGORITHM CALCULATIONS
CD2              = 1 CONVERGE WITHIN REQUIRED ACCURACY
CD2              > 1 ERROR ENCOUNTERED
CD2      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD2      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD2      N       NUMBER OF CONTROL VARIABLES
CD2      VLARGE  A VERY LARGE NUMBER
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      NVAR    (= N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      B       VARIABLE METRIC MATRIX. CORRESPONDS TO THE HESSIAN
CD5              MATRIX IN THE QUADRATIC FUNCTION
CD5      BDL     LOWER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD5              SEE EQN (1B)
CD5      BDU     UPPER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD5              SEE EQN (1B)
CD5      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD5              DIRECTION TIME THE STEP-LENGTH.
CD5      GM      NEGATIVE OF THE GRADIENT VECTOR (= -G). CORRESPONDS
CD5              TO B IN EQNS (4),(5),(6), AND (8).
CD5      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD5      MTOTAL  TOTAL NUMBER OF CONSTRAINTS (= M+2*NVAR). INCLUDES THE
CD5              INPUT CONSTRAINTS AND THE UPPER AND LOWER BOUNDARIES.
CD5      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD5              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN THE
CD5              BASIS
CD5     NVAR2    TWICE THE VALUE OF NVAR. (= 2*(N+1))
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD7
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION B      (MAXX  ,MAXX  )
C
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( B          ,  A      ( 1,1) )
      EQUIVALENCE ( MTOTAL     ,  KOMVMA (   4) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
      EQUIVALENCE ( NVAR2      ,  KOMVMA (   7) )
C
C***********************************************************************
C
C...  SET INITIAL VARIABLES AND ELEMENTS
C
      NVAR2 = NVAR*2
      MTOTAL = NVAR2 + M
      NC = MEQ + 1
C
      DO 10 I=1,N
        BDL(I) = -VLARGE
        BDU(I) = VLARGE
        DELTA(I) = ZERO
   10 CONTINUE
C
      BDL(NVAR) = ZERO
      DELTA(NVAR) = FP1
C
      IF (MEQ .GE. 1) THEN
        DO 20 I=1,MEQ
          IBASIS(I) = I + NVAR2
   20   CONTINUE
      ENDIF
C
      IBASIS(NC) = NVAR2
C
C...SET THE EXTRA ELEMENTS FOR GM AND B TO ALLOW FOR INFEASIBILITY
C
      GM(NVAR) = VLARGE
      IF (ISTAT .EQ. -1) THEN
        DO 30 I=1,NVAR
          B(I,NVAR) = ZERO
          B(NVAR,I) = ZERO
   30   CONTINUE
      ENDIF
C
      ISTAT = 0
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE STOVMU (G,M,N,X)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         STORES THE NECESSARY INFORMATION FOR THE GRADIENTS OF THE
CD1         LAGRANGE FUNCTION AND ALSO THE MULTIPLIERS FOR THE LINE
CD1         SEARCH FUNCTION
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      G       GRADIENT OF THE OBJECTIVE FUNCTION
CD2      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD2      N       NUMBER OF CONTROL VARIABLES
CD2      X       VECTOR OF CONTROL VARIABLES
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      GLF     GRADIENT OF THE LAGRANGE FUNCTION. SEE EQNS (3.2)
CD4              AND (3.3)
CD4      VLAM    VECTOR OF LAGRANGE MULTIPLIERS
CD5
CD5     COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      GDELTA  SCALAR PRODUCT OF GRADIENT OF THE OBJECTIVE FUNCTION
CD5              AND DELTA
CD5      GLFA    PREVIOUS GRADIENT OF THE LAGRANGE FUNCTION
CD5      OXA     SAVED VALUE OF THE VARIABLE XA
CD5      VMU     MULTIPLIERS FOR THE LINE SEARCH FUNCTION. EQN (4.2)
CD5      XA      INITIAL GUESS AT CONTROL VARIABLE VALUES
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      ABSLAM  ABSOLUTE VALUE OF LAGRANGE MULTIPLIERS (VLAM)
CD6      HALF    FLOATING POINT NUMBER 0.5
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXN    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS PLUS
CD6      NCOM    COMVMA ARRAY SIZE
CD6              THE ARTIFICIAL VARIABLE
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS CITED FROM "A FAST ALGORITHM FOR NONLINEARLY
CD7      CONSTRAINED OPTIMIZATION CALCULATIONS" BY M.J.D. POWELL
C
C***********************************************************************
C
      PARAMETER ( HALF   =  0.5 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXN   = MAXX - 1)
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION G      (MAXX  )
      DIMENSION X      (MAXX  )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAN / GLF    (MAXN  )
     1   ,              GLFA   (MAXN  )
     1   ,              OXA    (MAXN  )
     1   ,              XA     (MAXN  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
C
      EQUIVALENCE ( GDELTA     ,  COMVMA (   7) )
C
C***********************************************************************
C
C...FIRST THE SCALAR PRODUCT OF G AND DELTA IS CALCULATED. THE CONTROL
C...VALUE GUESSES AND THE GRADIENT ARE ALSO SAVED
C
      GDELTA = ZERO
C
      DO 10 I=1,N
        GDELTA = GDELTA + G(I)*DELTA(I)
        GLFA(I) = GLF(I)
        XA(I) = X(I)
        OXA(I) = XA(I)
   10 CONTINUE
C
      IF (M .GT. 0) THEN
C
C...    REVISE THE VMU VECTOR
C
        DO 20 I=1,M
          ABSLAM = ABS(VLAM(I))
          VMU(I) = ABIG(ABSLAM,HALF*(ABSLAM + VMU(I)))
   20   CONTINUE
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE UPDATB (N)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         UPDATES THE VARIABLE METRIC MATRIX B
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      N       NUMBER OF CONTROL VARIABLES
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH
CD4      GLF     GRADIENT OF THE LAGRANGE FUNCTION. SEE EQNS (3.2)
CD4              AND (3.3)
CD4      GLFA    PREVIOUS GRADIENT OF THE LAGRANGE FUNCTION
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      B       VARIABLE METRIC MATRIX. CORRESPONDS TO THE HESSIAN
CD5              MATRIX IN THE QUADRATIC FUNCTION
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      BDDBD   BDELTA DIVIDED BY DBD
CD6      BDELTA  PRODUCT OF THE VARIABLE METRIC MATRIX B AND DELTA
CD6      BFACTR  B MATRIX FACTOR. CHOSEN EMPIRICALLY TO USE IN THE
CD6              B-F-G-S FORMULA. SEE EQN (3.6)
CD6      DBD     SCALAR PRODUCT OF DELTA*B*DELTA. SEE EQN (3.6)
CD6      DBDFAC  PRODUCT OF DBD AND BFACTR. SEE EQN (3.6)
CD6      DELETA  SCALAR PRODUCT OF DELTA AND ETA. SEE EQN (3.8)
CD6      DELGAM  SCALAR PRODUCT OF DELTA AND GAMMA WHERE GAMMA
CD6              CORRESPONDS TO THE DIFFERENCE IN GRADIENTS. EQN (3.3)
CD6      DIFGRD  DIFFERENCE BETWEEN THE CURRENT AND PREVIOUS GRADIENT OF
CD6              THE LAGRANGE FUNCTION. CORRESPONDS TO GAMMA IN EQN(3.3)
CD6      ETA     REPLACES DIFGRD (I.E. GAMMA) IN THE B-F-G-S FORMULA.
CD6              SEE EQN (3.5)
CD6      ETADE   ETA DIVIDED BY DELETA
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      MAXN    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6              AND THE ARTIFICIAL VARIABLE
CD6      THETA   THE THETA OF EQN (3.6) AND (3.7)
CD6      THCOMP  COMPLIMENT OF THETA (I.E. 1-THETA)
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS CITED FROM "A FAST ALGORITHM FOR NONLINEARLY
CD7      CONSTRAINED OPTIMIZATION CALCULATIONS" BY M.J.D. POWELL
C
C***********************************************************************
C
      PARAMETER ( BFACTR =  0.2 )
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXN   = MAXX - 1)
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION B      (MAXX  ,MAXX  )
      DIMENSION BDELTA (MAXN  )
      DIMENSION DIFGRD (MAXN  )
      DIMENSION ETA    (MAXN  )
C
      COMMON   /CVMAN / GLF    (MAXN  )
     1   ,              GLFA   (MAXN  )
     1   ,              OXA    (MAXN  )
     1   ,              XA     (MAXN  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
C
      EQUIVALENCE ( B          ,  A      ( 1,1) )
C
C***********************************************************************
C
C...CALCULATE THE NECESSARY MATRICES AND INFORMATION TO REVISE THE
C...B MATRIX
C
      DELGAM = ZERO
      DBD = ZERO
C
      DO 20 I=1,N
        DIFGRD(I) = GLF(I)-GLFA(I)
        BDELTA(I) = ZERO
        DO 10 J=1,N
          BDELTA(I) = BDELTA(I) + B(I,J)*DELTA(J)
   10   CONTINUE
C
        DELGAM = DELGAM + DELTA(I)*DIFGRD(I)
        DBD = DBD + DELTA(I)*BDELTA(I)
   20 CONTINUE
C
C...  CALCULATE THE VECTOR ETA FOR THE B-F-G-S FORMULA
C
      DBDFAC = BFACTR*DBD
C
      IF (DELGAM .LT. DBDFAC) THEN
        THETA = (DBD - DBDFAC)/(DBD-DELGAM)
        THCOMP = FP1-THETA
        DO 30 I=1,N
          ETA(I) = THETA*DIFGRD(I) + THCOMP*BDELTA(I)
   30   CONTINUE
C
        DELETA = DBDFAC
      ELSE
        DO 40 I=1,N
          ETA(I) = DIFGRD(I)
   40   CONTINUE
        DELETA = DELGAM
      ENDIF
C
C...REVISE THE B MATRIX AS PER EQN (3.8)
C
      DO 60 I=1,N
        BDDBD = BDELTA(I)/DBD
        ETADE = ETA(I)/DELETA
C
        DO 50 J=I,N
          B(I,J) = B(I,J) - BDDBD*BDELTA(J) + ETADE*ETA(J)
          B(J,I) = B(I,J)
   50   CONTINUE
   60 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE VERTEX (MEQ)
CD0
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         OBTAINS AN INITIAL FEASIBLE VERTEX (A POINT FORMED BY THE
CD1         INTERSECTION OF NVAR INEQUALITIES) AND WORKS BY EXCHANGING
CD1         CONSTRAINTS IN A TRIAL VERTEX
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD2      NVIOL   NUMBER OF VIOLATED CONSTRAINTS
CD2      SNORM   SUM OF THE VIOLATED CONSTRAINT NORMALS
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD3      ICROW   INDEX OF THE ROW OF CINVRS WHICH IS TO BE REMOVED FROM
CD3              THE BASIS
CD3      LSIZE   TEMPORARY ARRAY SIZE OF A VECTOR BEING MULTIPLIED
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      BDL     LOWER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1A)
CD4      BDU     UPPER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE
CD4              SEE EQN (1A)
CD4      CINVRS  INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD4              CONSTRAINTS IN THE BASIS. SEE EQN (2)
CD4      D       RIGHT HAND SIDE OF CONSTRAINT EQUATION. SEE EQN (2)
CD4      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD4      MTOTAL  TOTAL NUMBER OF CONSTRAINTS (= M+2*NVAR). INCLUDES THE
CD4              INPUT CONSTRAINTS AND THE UPPER AND LOWER BOUNDARIES
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN THE
CD4              BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD4      NVAR2   TWICE THE VALUE OF NVAR. (= 2*(N+1))
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD5      IACTIV  STATUS OF EACH OF THE MTOTAL CONSTRAINTS IN THE BASIS
CD5              =-1 EQUALITY
CD5              = 0 ACTIVE
CD5              = 1 INACTIVE
CD5              = 2 VIOLATED
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      MAXACT  MAXIMUM SIZE OF ARRAY DEALING WITH IACTIV ARRAY
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS PLUS
CD6      NKOM    KOMVMA ARRAY SIZE
CD6              THE ARTIFICIAL VARIABLE
CD6      RHS     RIGHT HAND SIDE OF THE CONSTAINT EQNS IN THE BASIS
CD6      UICJ    PRODUCT OF THE ITH ROW OF CINVRS AND THE NORMAL VECTORS
CD6              OF CONSTRAINTS IN THE BASIS
CD6      UICMAX  MAXIMUM VALUE OF UICJ
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "THE CALCULATION OF FEASIBLE
CD7      POINTS FOR LINEARLY CONSTRAINED OPTIMIZATION PROBLEMS" AERE-
CD7      R.6354 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXACT = MAXX*2 + MAXC )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION RHS    (MAXX  )
      DIMENSION SNORM  (MAXC  )
C
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMACT/ IACTIV (MAXACT)
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( MTOTAL     ,  KOMVMA (   4) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
      EQUIVALENCE ( NVAR2      ,  KOMVMA (   7) )
      EQUIVALENCE ( ICROW      ,  KOMVMA (  10) )
C
C***********************************************************************
C
C...INITIALIZE THE STATUS OF THE CONSTRAINTS AS INACTIVE
C
      DO 10 I=1,MTOTAL
        IACTIV(I) = 1
   10 CONTINUE
C
C...THE INVERSE MATRIX CORRESPONDING TO THE CONSTRAINTS IN THE BASIS
C...ARE CALCULATED
C
      IF (NC .EQ. 0) THEN
        CALL NODESG
      ELSE
        CALL DESIG (MEQ)
      ENDIF
C
C...THE RIGHT HAND SIDE OF THE CONSTRAINTS IN THE BASIS ARE CALCULATED
C
      DO 20 I=1,NVAR
        IF (IBASIS(I) .GT. NVAR2) THEN
C
C...      EQUALITY OR INEQUALITY CONSTRAINT IN THE BASIS
C
          RHS(I) = D(IBASIS(I)-NVAR2)
        ELSEIF (IBASIS(I) .GT. NVAR) THEN
C
C...      UPPER BOUNDARY CONSTRAINT IN THE BASIS
C
          RHS(I) = -BDU(IBASIS(I)-NVAR)
        ELSE
C
C...      LOWER BOUNDARY CONSTRAINT IN THE BASIS
C
          RHS(I) = BDL(IBASIS(I))
        ENDIF
   20 CONTINUE
C
C...CALCULATE THE POSITION OF THE VERTEX
C
      DO 30 I=1,NVAR
        DELTA(I) = DOTPRD(NVAR,CINVRS(1,I),NVAR,1,RHS(1),NVAR,1)
   30 CONTINUE
C
C...A TRIAL VERTEX HAS BEEN FOUND. IT MUST BE DETERMINED IF IT IS A
C...FEASIBLE POINT
C
   40 CONTINUE
      CALL CKVRTX (NVIOL,SNORM)
C
      IF (NVIOL .NE. 0) THEN
C
C...    POSSIBLE DIRECTIONS OF SEARCH ARE ROWS OF CINVRS OBTAINED BY
C...    REMOVING A CONSTRAINT. CALCULATE THE OPTIMUM DIRECTION.
C
        UICMAX = ZERO
        LSIZE = 1 + MAXX*(NVAR - 1)
        DO 50 I=1,NVAR
          IF (IACTIV(IBASIS(I)) .NE. -1) THEN
            UICJ = DOTPRD (NVAR,CINVRS(I,1),LSIZE,MAXX,SNORM(1),NVAR,1)
            IF (UICJ .GT. UICMAX) THEN
              UICMAX = UICJ
              ICROW = I
            ENDIF
          ENDIF
   50   CONTINUE
C
        IF (UICMAX .LE. ZERO) THEN
          CALL VMAMSG(11,0,ZERO)
        ELSE
          CALL NEWBAS (ICROW)
          GO TO 40
        ENDIF
      ENDIF
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE VMACO (N,M,MEQ,X,F,G,C,CNORML,NSIZE,MAXFUN,
     1                  ACC,ISTAT,IPRINT)
CD0
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         VARIABLE METRIC ALGORITHM FOR CONSTRAINED OPTIMIZATION
CD1
CD1         CALCULATES THE LEAST VALUE OF A FUNCTION OF SEVERAL VARIABLES
CD1         SUBJECT TO EQUALITY AND INEQUALITY CONSTRAINTS.
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      ACC     CONTROLS THE FINAL ACCURACY FOR CONVERGENCE. CON-
CD2              VERGENCE OCCURS WHEN THE VALUE FOR THE CHANGE IN
CD2              OBJECTIVE FUNCTION PLUS SUITABLY WEIGHTED MULTI-
CD2              PLIERS OF THE CONSTRAINT FUNCTION DIFFERS FROM
CD2              ITS PREVIOUS VALUE BY AT MOST ACC.
CD2      C       VECTOR OF CONSTRAINT VARIABLES.
CD2      CNORML  MATRIX OF CONSTRAINT NORMALS
CD2      F       VALUE OF THE OBJECTIVE FUNCTION (OPTIMIZED VARIABLE)
CD2      G       GRADIENT OF THE OBJECTIVE FUNCTION
CD2      IPRINT  FLAG TO INDICATE AMOUNT OF PRINTOUT
CD2              = 0 NO PRINTOUT
CD2              = 1 VALUES OF X, C, AND F ARE PRINTED
CD2              = 2 DEBUGGING PRINTOUT
CD2      ISTAT   INDEX TO INDICATE THE STATUS OF THE OPTIMIZATION
CD2              ALGORITHM.
CD2              =-2 INITIAL PASS THROUGH VMACO WITH THE B MATRIX
CD2                  ALREADY INITIALIZED
CD2              =-1 INITIAL PASS THROUGH VMACO
CD2              = 0 CONTINUE ALGORITHM CALCULATIONS
CD2              = 1 CONVERGE WITHIN REQUIRED ACURRACY
CD2              > 1 ERROR ENCOUNTERED
CD2      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD2      MAXFUN  MAXIMUM NUMBER OF TIMES VMACO CAN BE CALLED
CD2      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD2      N       NUMBER OF CONTROL VARIABLES
CD2      NSIZE   CORRESPONDS TO THE MAXIMUM DIMENSION OF THE INNERMOST
CD2              ARRAY SIZE FROM THE CNORML MATRIX (I.E. MAX CONTROL
CD2              VARIABLE SIZE). THIS ALLOWS CORRECT COMMUNICATION
CD2              BETWEEN THE CNORML AND CN MATRIX.
CD2      X       VECTOR OF CONTROL VARIABLES
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      BFLAG   FLAG FOR WHETHER THE B MATRIX SHOULD BE UPDATED
CD3      CNORML  MATRIX OF CONSTRAINT NORMALS
CD3      ISTAT   INDEX TO INDICATE THE STATUS OF THE OPTIMIZATION
CD3              ALGORITHM
CD3              =-1 INITIAL PASS THROUGH VMACO
CD3              = 0 CONTINUE ALGORITHM CALCULATIONS
CD3              = 1 CONVERGE WITHIN REQUIRED ACCURACY
CD3              > 1 ERROR ENCOUNTERED
CD3      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD3      MAXFUN  MAXIMUM NUMBER OF TIMES VMACO CAN BE CALLED
CD3      N       NUMBER OF CONTROL VARIABLES
CD3      NVMITR  VALUE OF NVMA AT THE START OF AN ITERATION
CD3      X       VECTOR OF CONTROL VARIABLES
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      CN      MATRIX OF CONSTRAINT NORMALS
CD4      GDELTA  SCALAR PRODUCT OF GRADIENT OF THE OBJECTIVE FUNCTION
CD4              AND DELTA
CD4      NVMA    NUMBER OF CALLS TO VMACO
CD4      VLAM    VECTOR OF LAGRANGE MULTIPLIERS
CD4      X       VECTOR OF CONTROLS PLUS AND ARTIFICIAL VARIABLE
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      ITER    COUNTER FOR NUMBER OF ITERATIONS COMPLETED
CD5      IBUGFG  FLAG TO INDICATE EXTRA PRINTOUT FOR DEBUGGING PURPOSES
CD5              (= 0) NO EXTRA PRINTOUT
CD5              (= 1) EXTRA PRINTOUT
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FDELTA  ABSOLUTE CHANGE IN THE OBJECTIVE FUNCTION F. CHECKED
CD6              AGAINST ACCURACY DESIRED FOR CONVERGENCE CRITERIA
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXCM1  VALUE OF MAXC SUBTRACT 1
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH THE CONTROLS PLUS
CD6              THE ARTIFICIAL VARIABLE
CD6      MAXXM1  VALUE OF MAXX SUBTRACT 1
CD6      NCOM    COMVMA ARRAY SIZE
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      NVMLIN  CONTROLS THE ERROR RETURN FROM THE LINE SEARCH
CD6              OBJECTIVE FUNCTION. SEE EQN (4.2)
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      ALGORITHM AND EQUATIONS CITED BASED ON THE PAPER "A FAST
CD7      ALGORITHM FOR NONLINEARLY CONSTRAINED OPTIMIZATION
CD7      CALCULATIONS" BY M.J.D POWELL PRESENTED AT THE 1977 DUNDEE
CD7      CONFERENCE ON NUMERICAL ANALYSIS
C
C***********************************************************************
C
      LOGICAL BFLAG
C
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( MAXCM1 = MAXC-1 )
      PARAMETER ( MAXXM1 = MAXX-1 )
      PARAMETER ( NCOM   =   15 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( NVMLIN =    5 )
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION C      (MAXC  )
      DIMENSION CNORML (NSIZE ,MAXC  )
      DIMENSION G      (MAXX  )
      DIMENSION X      (MAXX  )
C
      COMMON   /COMVMA/ COMVMA (NCOM  )
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( GDELTA     ,  COMVMA (   7) )
      EQUIVALENCE ( FDELTA     ,  COMVMA (  10) )
      EQUIVALENCE ( ITER       ,  KOMVMA (   2) )
      EQUIVALENCE ( NVMA       ,  KOMVMA (   8) )
      EQUIVALENCE ( NVMITR     ,  KOMVMA (   9) )
      EQUIVALENCE ( IBUGFG     ,  KOMVMA (  11) )
C
C***********************************************************************
C
C...CHECK N, M, AND MEQ ARE LOGICAL VALUES
C
      IF (N .LE. 0) THEN
        ISTAT=2
        CALL VMAMSG (ISTAT,0,ZERO)
        GO TO 999
      ELSEIF (M .LT. MEQ) THEN
        ISTAT=3
        CALL VMAMSG (ISTAT,0,ZERO)
        GO TO 999
      ELSEIF (MEQ .LT. 0) THEN
        ISTAT = 4
        CALL VMAMSG (ISTAT,0,ZERO)
        GO TO 999
      ELSEIF (N .GT. MAXXM1) THEN
        ISTAT = 12
        CALL VMAMSG (ISTAT, MAXXM1, ZERO)
        GO TO 999
      ELSEIF (M .GT. MAXCM1) THEN
        ISTAT = 13
        CALL VMAMSG (ISTAT, MAXCM1, ZERO)
        GO TO 999
      ENDIF
C
C...EQUATE CNORML TO CN TO ALLOW FOR ADDING AN EXTRA
C...VARIABLE FOR INFEASIBILITY
C
      IF (M .GE. 1) THEN
        DO 5 I=1,N
          DO 4 J=1,M
            CN(I,J)= CNORML(I,J)
    4     CONTINUE
    5   CONTINUE
      ENDIF
C
      IF (ISTAT.LE.-1) THEN
C
C...    CURRENTLY ON THE INITIAL ITERATION PASS
C
        CALL INITVM (ISTAT,N,M)
        GO TO 20
      ENDIF
C
      IF (IPRINT .EQ. 2) THEN
        IF (IBUGFG .EQ. 1) THEN
          CALL VMAOUT (C,F,G,M,N,X,IPRINT)
          IBUGFG = 0
        ENDIF
      ENDIF
C
      IF (NVMA .GE. MAXFUN) THEN
        ISTAT = 5
        CALL VMAMSG (ISTAT,MAXFUN,ZERO)
        GO TO 999
      ELSEIF (NVMA .GT. NVMITR+NVMLIN) THEN
        ISTAT = 6
        CALL VMAMSG (ISTAT,NVMLIN,ZERO)
        GO TO 999
      ENDIF
C
C...CALCULATE THE LINE SEARCH OBJECTIVE FUNCTION
C
   10 CONTINUE
      BFLAG = .FALSE.
C
      CALL LNSRCH (C,F,G,ISTAT,M,MAXFUN,MEQ,N,NVMITR,NVMLIN,
     1             X,BFLAG,IPRINT)
C
      IF (BFLAG .AND. ISTAT .LE. 0) THEN
C
C...    UPDATE THE B MATRIX
C
        CALL UPDATB (N)
      ELSE
        IBUGFG = 1
        GO TO 999
      ENDIF
C
   20 CONTINUE
      ITER = ITER + 1
C
      IF (IPRINT .GE. 1) THEN
        CALL VMAOUT (C,F,G,M,N,X,IPRINT)
        IF (IPRINT .EQ. 2) THEN
          IBUG = 2
          CALL BUGAID (IBUG,N,XDUMMY)
        ENDIF
      ENDIF
C
      CALL VMASET (C,G,ISTAT,M,MEQ,N,IPRINT)
C
      IF (ISTAT .GE. 2) THEN
        GO TO 999
      ENDIF
C
      NVMITR = NVMA
      CALL GRADLF (G,M,N,IPRINT)
      CALL STOVMU (G,M,N,X)
C
C...CALCULATE WHETHER THE ITERATION HAS CONVERGED ON THE OPTIMAL SOLUTION
C
      FDELTA = ABS(GDELTA)
      DO 40 I=1,M
        FDELTA = FDELTA + ABS(VLAM(I)*C(I))
   40 CONTINUE
C
      IF (FDELTA .LE. ACC) THEN
C
C...    CONVERGENCE HAS OCURRED
C
        ISTAT = 1
        IF (IPRINT .GE. 1) THEN
          CALL VMAMSG (ISTAT,0,ZERO)
        ENDIF
      ELSE
C
C...    DO A LINE SEARCH FOR NEW VALUES
C
        GO TO 10
      ENDIF
C
  999 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
      SUBROUTINE VMAMSG (ICODE,NUM,XNUM)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         OUTPUTS DIAGNOSTIC MESSAGES
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      ICODE   MESSAGE CODE NUMBER
CD2      NUM     INTEGER NUMBER
CD2      XNUN    REAL NUMBER
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      NONE
CD4
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      NONE
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      BORDER  DIVIDING LINE OF MESSAGES
CD6      IECHO   FILE UNIT NUMBER TO WHICH ERROR MESSAGES ARE WRITTEN
CD6      MESSGE  ARRAY OF MESSAGES
CD6      NEXTRA  NUMBER OF EXTRA 72 WORD CHARACTER MESSAGES
CD6      NUMMSG  SIZE OF ARRAY MESSAGE
CD6
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      NONE
C
C***********************************************************************
C
      CHARACTER BORDER * 72
      CHARACTER MESSGE * 72
      CHARACTER DESIG  *  5
      CHARACTER RESET  *  5
C
      PARAMETER ( IECHO  =    6 )
      PARAMETER ( NUMMSG =   15 )
      PARAMETER ( NEXTRA =    2 )
C
      DIMENSION MESSGE (NUMMSG,NEXTRA)
C
      DATA BORDER       ( 1:36) /'------------------------------------'/
      DATA BORDER       (37:72) /'------------------------------------'/
C
      DATA DESIG        ( 1:5 ) /'DESIG'/
      DATA RESET        ( 1:5 ) /'RESET'/
C
      DATA MESSGE( 1, 1)( 1:36) /'THE VARIABLE METRIC ALGORITHM HAS DE'/
      DATA MESSGE( 1, 1)(37:72) /'TERMINED AN OPTIMAL SOLUTION WITHIN '/
      DATA MESSGE( 1, 2)( 1:36) /'THE ACCURACY DESIRED                '/
      DATA MESSGE( 1, 2)(37:72) /'                                    '/
      DATA MESSGE( 2, 1)( 1:36) /'*** ERROR *** THE INPUT NUMBER OF CO'/
      DATA MESSGE( 2, 1)(37:72) /'NTROL VARIABLES IS LESS THAN ONE    '/
      DATA MESSGE( 2, 2)( 1:36) /'                                    '/
      DATA MESSGE( 2, 2)(37:72) /'                                    '/
      DATA MESSGE( 3, 1)( 1:36) /'*** ERROR *** THE INPUT NUMBER OF CO'/
      DATA MESSGE( 3, 1)(37:72) /'NSTRAINTS IS LESS THAN THE NUMBER OF'/
      DATA MESSGE( 3, 2)( 1:36) /' EQUALITY CONSTRAINTS               '/
      DATA MESSGE( 3, 2)(37:72) /'                                    '/
      DATA MESSGE( 4, 1)( 1:36) /'*** ERROR *** THE INPUT NUMBER OF EQ'/
      DATA MESSGE( 4, 1)(37:72) /'UALITY CONSTRAINTS IS LESS THEN ZERO'/
      DATA MESSGE( 4, 2)( 1:36) /'                                    '/
      DATA MESSGE( 4, 2)(37:72) /'                                    '/
      DATA MESSGE( 5, 1)( 1:36) /'*** WARNING *** MAXIMUM NUMBER OF CA'/
      DATA MESSGE( 5, 1)(37:72) /'LLS (=XXXX) TO VMACO HAS BEEN       '/
      DATA MESSGE( 5, 2)( 1:36) /'EXCEEDED                            '/
      DATA MESSGE( 5, 2)(37:72) /'                                    '/
      DATA MESSGE( 6, 1)( 1:36) /'*** ERROR *** LINE SEARCH REQUIRES X'/
      DATA MESSGE( 6, 1)(37:72) /'XXX CALLS TO VMACO                  '/
      DATA MESSGE( 6, 2)( 1:36) /'                                    '/
      DATA MESSGE( 6, 2)(37:72) /'                                    '/
      DATA MESSGE( 7, 1)( 1:36) /'*** ERROR *** AN UPHILL SEARCH DIREC'/
      DATA MESSGE( 7, 1)(37:72) /'TION HAS BEEN ENCOUNTERED. USUALLY  '/
      DATA MESSGE( 7, 2)( 1:36) /'THIS IS DUE TO LOSS OF ACCURACY IN T'/
      DATA MESSGE( 7, 2)(37:72) /'HE QUADRATIC PROGRAMMING CALCULATION'/
C
      DATA MESSGE( 8, 1)( 1:36) /'*** ERROR *** THE GIVEN CONSTRAINT S'/
      DATA MESSGE( 8, 1)(37:72) /'EEMS TO BE INCONSISTENT             '/
      DATA MESSGE( 8, 2)( 1:36) /'                                    '/
      DATA MESSGE( 8, 2)(37:72) /'                                    '/
      DATA MESSGE( 9, 1)( 1:36) /'*** ERROR *** AN ARTIFICIAL BOUND IS'/
      DATA MESSGE( 9, 1)(37:72) /' ACTIVE. THE PREDICTED CHANGE IN THE'/
      DATA MESSGE( 9, 2)( 1:36) /'VARIABLE EXCEEDS XXXXXX             '/
      DATA MESSGE( 9, 2)(37:72) /'                                    '/
      DATA MESSGE(10, 1)( 1:36) /'*** WARNING *** THE SOLUTION MAY DEG'/
      DATA MESSGE(10, 1)(37:72) /'ENERATE INTO A LOCAL MINIMUM        '/
      DATA MESSGE(10, 2)( 1:36) /'                                    '/
      DATA MESSGE(10, 2)(37:72) /'                                    '/
      DATA MESSGE(11, 1)( 1:36) /'*** WARNING *** NO FEASIBLE VERTEX F'/
      DATA MESSGE(11, 1)(37:72) /'OUND                                '/
      DATA MESSGE(11, 2)( 1:36) /'                                    '/
      DATA MESSGE(11, 2)(37:72) /'                                    '/
      DATA MESSGE(12, 1)( 1:36) /'*** ERROR *** THE MAXIMUM SIZE OF TH'/
      DATA MESSGE(12, 1)(37:72) /'E CONTROL ARRAY (NSIZE) IS GREATER  '/
      DATA MESSGE(12, 2)( 1:36) /'THAN XXXX                           '/
      DATA MESSGE(12, 2)(37:72) /'                                    '/
      DATA MESSGE(13, 1)( 1:36) /'*** ERROR *** THE MAXIMUM SIZE OF TH'/
      DATA MESSGE(13, 1)(37:72) /'E CONSTRAINT ARRAY ( M ) IS         '/
      DATA MESSGE(13, 2)( 1:36) /'GREATER THAN XXXX                   '/
      DATA MESSGE(13, 2)(37:72) /'                                    '/
      DATA MESSGE(14, 1)( 1:36) /'*** ERROR *** QUADRATIC SEARCH REQUI'/
      DATA MESSGE(14, 1)(37:72) /'RES MORE THAN XXXX LOOPS.           '/
      DATA MESSGE(14, 2)( 1:36) /'LIMIT OR ADJUST THE APPROPRIATE CONT'/
      DATA MESSGE(14, 2)(37:72) /'ROLS                                '/
      DATA MESSGE(15, 1)( 1:36) /'*** WARNING *** RANK OF AN INVERTED '/
      DATA MESSGE(15, 1)(37:72) /'MATRIX IN SUBROUTINE XXXXX HAS      '/
      DATA MESSGE(15, 2)( 1:36) /'DECREASED                           '/
      DATA MESSGE(15, 2)(37:72) /'                                    '/
C
C***********************************************************************
C
C...SET PROPER PARAMETER VALUES
C
      IF (ICODE .EQ. 5) THEN
        WRITE (MESSGE(5,1)(43:46), 10) NUM
      ELSEIF (ICODE .EQ. 6) THEN
        WRITE (MESSGE(6,1)(36:39), 10) NUM
      ELSEIF (ICODE .EQ. 9) THEN
        WRITE (MESSGE(9,2)(18:33), 20) XNUM
      ELSEIF (ICODE .EQ. 12 ) THEN
        WRITE (MESSGE(12,2)(6:9), 10) NUM
      ELSEIF (ICODE .EQ. 13 ) THEN
        WRITE (MESSGE(13,2)(14:17) , 10) NUM
      ELSEIF (ICODE .EQ. 14 ) THEN
        WRITE (MESSGE(14,1)(51:54) , 10) NUM
      ELSEIF (ICODE .EQ. 15 ) THEN
        IF (NUM .EQ. 1) THEN
            MESSGE(15,1)(58:62) = DESIG
        ELSEIF (NUM .EQ. 2) THEN
            MESSGE(15,1)(58:62) = RESET
        ENDIF
      ENDIF
C
C...WRITE DIAGNOSTIC MESSAGE
C
      WRITE (IECHO, 30) BORDER,MESSGE(ICODE,1),MESSGE(ICODE,2),ICODE,
     1                  BORDER
C
C***********************************************************************
C
      RETURN
C
   10 FORMAT (I4)
   20 FORMAT (1PE16.7)
   30 FORMAT (' ',/,1X,A72,//,1X,A72,/,1X,A72,
     1        /,1X,'*** ISTAT = ',I2,' ***',/,1X,A72,/)
C
      END
      SUBROUTINE VMAOUT (C,F,G,M,N,X,IPRINT)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         PROVIDES THE PRINTOUT OF THE OBJECTIVE FUNCTION, CONTROL
CD1         VARIABLE, AND CONSTRAINT VARIABLE VALUES UPON EACH RETURN
CD1         FROM VMACO
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      C       VECTOR OF CONSTRAINT VARIABLES
CD2      G       GRADIENT OF THE OBJECTIVE FUNCTION
CD2      IPRINT  FLAG TO INDICATE AMOUNT OF PRINTOUT
CD2              =0 NO PRINTOUT
CD2              =1 VALUES OF X, C, AND F ARE PRINTED
CD2              =2 DEBUGGING PRINTOUT
CD2      F       VALUE OF THE OBJECTIVE FUNCTION (OPTIMIZED VARIABLE)
CD2      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD2      N       NUMBER OF CONTROL VARIABLES
CD2      X       VECTOR OF CONTROL VARIABLES
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      NONE
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      ITER    COUNTER FOR NUMBER OF ITERATIONS COMPLETED
CD4      NVMA    NUMBER OF CALLS TO VMACO
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      NONE
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      IECHO   FILE UNIT NUMBER TO WHICH PRINTOUT IS WRITTEN
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6      NKOM    KOMVMA ARRAY SIZE
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      NONE
C
C***********************************************************************
C
      PARAMETER ( IECHO  =    6 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NKOM   =   15 )
C
      DIMENSION C      (MAXC  )
      DIMENSION G      (MAXX  )
      DIMENSION X      (MAXX  )
C
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
C
      EQUIVALENCE ( ITER       ,  KOMVMA (   2) )
      EQUIVALENCE ( NVMA       ,  KOMVMA (   8) )
C
C***********************************************************************
C
      WRITE (IECHO,  5)
      WRITE (IECHO, 10) ITER, NVMA
      WRITE (IECHO, 20)
      CALL COLUM3 (6,N,X,NDUMMY)
C
      IF (M .GE. 1) THEN
        WRITE (IECHO, 40)
        CALL COLUM3 (7,M,C,NDUMMY)
      ENDIF
C
      WRITE (IECHO, 60) F
      IF (IPRINT .EQ. 2) THEN
        WRITE (IECHO,70)
        CALL COLUM3 (8,N,G,NDUMMY)
        WRITE (IECHO,80)
        DO 3 I=1,M
          CALL COLUM3 (9,N,XDUMMY,I)
    3   CONTINUE
      ENDIF
      WRITE (IECHO,  5)
C
C***********************************************************************
C
      RETURN
C
    5 FORMAT (1X,72(1H*))
   10 FORMAT (/5X, 12HITERATIONS =, I5, 5X,
     1        26HNUMBER OF CALLS TO VMACO =, I5)
   20 FORMAT (/1X,42HTHE VECTOR OF CONTROL VARIABLE VALUES ARE:/)
   40 FORMAT (/1X, 36HTHE VECTOR OF CONSTRAINT VALUES ARE:/)
   60 FORMAT (/1X, 32HTHE OBJECTIVE FUNCTION VALUE IS://5X, 3HF= ,
     1        1PE16.7)
   70 FORMAT (/2X,34HTHE VECTOR OF GRADIENT VALUES ARE:/)
   80 FORMAT (/2X,23HTHE JACOBIAN MATRIX IS:/)
C
      END
      SUBROUTINE VMASET (C,G,ISTAT,M,MEQ,N,IPRINT)
CDO
CD0   IDENTIFICATION :
CD0
CD0         PROGRAMMER -  J.D. FRICK, MDAC/HOUSTON
CD0         VERSION    -  JAN 1986
CD1
CD1   PURPOSE :
CD1
CD1         SETS UP AND CALCULATES VARIABLES NEEDED BY VMACO AND THE
CD1         QUADRATIC PROGRAMMING ROUTINES
CD2
CD2   CALLING ARGUMENT INPUT :
CD2
CD2      NAME    DEFINITION
CD2
CD2      C       VECTOR OF CONSTRAINT VARIABLES
CD2      G       GRADIENT OF THE OBJECTIVE FUNCTION
CD2      IPRINT  FLAG TO INDICATE AMOUNT OF PRINTOUT
CD2              = 0 NO PRINTOUT
CD2              = 1 VALUES OF X, C, AND F ARE PRINTED
CD2              = 2 DEBUGGING PRINTOUT
CD2      ISTAT   INDEX TO INDICATE THE STATUS OF THE OPTIMIZATION
CD2              ALGORITHM
CD2              =-2 INITIAL PASS THROUGH VMACO WITH THE B MATRIX
CD2                  ALREADY INITIALIZED
CD2              =-1 INITIAL PASS THROUGHT VMACO
CD2              = 0 CONTINUE ALGORITHM CALCULATIONS
CD2              = 1 CONVERGE WITHIN REQUIRED ACCURACY
CD2              > 1 ERROR ENCOUNTERED
CD2      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD2      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD2      N       NUMBER OF CONTROL VARIABLES
CD3
CD3   CALLING ARGUMENT OUTPUT :
CD3
CD3      NAME    DEFINITION
CD3
CD3      ISTAT   INDEX TO INDICATE THE STATUS OF THE OPTIMIZATION
CD3              ALGORITHM
CD2              =-2 INITIAL PASS THROUGH VMACO WITH THE B MATRIX
CD2                  ALREADY INITIALIZED
CD3              =-1 INITIAL PASS THROUGHT VMACO
CD3              = 0 CONTINUE ALGORITHM CALCULATIONS
CD3              = 1 CONVERGE WITHIN REQUIRED ACCURACY
CD3              > 1 ERROR ENCOUNTERED
CD3      M       TOTAL NUMBER OF CONSTRAINTS (I.E. CONSTRAINT EQUATIONS)
CD3      MEQ     TOTAL NUMBER OF EQUALITY CONSTRAINTS
CD3      N       NUMBER OF CONTROL VARIABLES
CD3      VLARGE  A VERY LARGE NUMBER
CD4
CD4   COMMON VARIABLE INPUT
CD4
CD4      NAME    DEFINITION
CD4
CD4      CSTAR   INVERSE MATRIX CORRESPONDING TO THE NORMALS OF THE
CD4              CONSTRAINTS IN THE BASIS. SEE EQN (11)
CD4              EQUIVALENT TO CINVRS IN THE ROUTINE VERTEX
CD4      DELTA   CHANGE IN THE CONTROL VARIABLE. CORRESPONDS TO SEARCH
CD4              DIRECTION TIMES THE STEP-LENGTH.
CD4      GHAT    CURRENT GRADIENT VECTOR OF THE OBJECTIVE FUNCTION.
CD4              SEE EQN (9)
CD4      IBASIS  INDEX NUMBER OF THE CONSTRAINTS IN THE BASIS
CD4      NC      NUMBER OF CONSTRAINTS IN THE BASIS. INITIALLY
CD4              REPRESENTS THE NUMBER OF DESIGNATED CONSTRAINTS IN
CD4              THE BASIS
CD4      NVAR    (=N+1) NUMBER OF CONTROLS PLUS THE ARTIFICIAL VARIABLE
CD4      NVAR2   TWICE THE VALUE OF NVAR. ( = 2*(N+1) )
CD4      VLAM    VECTOR OF LAGRANGE MULTIPLIERS
CD5
CD5   COMMON VARIABLE OUTPUT
CD5
CD5      NAME    DEFINITION
CD5
CD5      BDU     UPPER BOUNDS TO THE CONTROLS AND ARTIFICIAL VARIABLE.
CD5              SEE EQN (1B)
CD5      CN      MATRIX OF CONSTRAINT NORMALS
CD5      D       RIGHT HAND SIDE OF CONSTRAINT EQUATION. SEE EQN (1C)
CD5      GM      NEGATIVE OF THE GRADIENT VECTOR (= -G). CORRESPONDS
CD5              TO B IN EQNS (4),(5),(6), AND (8).
CD5      VLAM    VECTOR OF LAGRANGE MULTIPLIERS
CD6
CD6   INTERNAL VARIABLES
CD6
CD6      NAME    DEFINITION
CD6
CD6      FEASP   SCALING FACTOR USED TO ACHIEVE FEASIBILITY. CORRESPONDS
CD6              TO GREEK LETTER XI, EQN (2.7) IN POWELL'S PAPER.
CD6      FP1     FLOATING POINT NUMBER 1.0
CD6      IFEAS   FLAG USED TO MARK FEASIBILITY CONDITIONS
CD6      MAXC    MAXIMUM SIZE OF ARRAYS DEALING WITH CONSTRAINT VARIABLE
CD6      MAXX    MAXIMUM SIZE OF ARRAYS DEALING WITH CONTROL VARIABLES
CD6      NKOM    KOMVMA ARRAY SIZE
CD6      VSMALL  A VERY SMALL POSITIVE NUMBER
CD6      ZERO    FLOATING POINT NUMBER 0.0
CD7
CD7   LIMITATIONS/ASSUMPTIONS
CD7
CD7      EQUATIONS AND ALGORITHM CITED FROM "A GENERAL QUADRATIC
CD7      ALGORITHM" AERE-TP.401 BY R. FLETCHER
C
C***********************************************************************
C
      PARAMETER ( FEASP  =  0.9 )
      PARAMETER ( FP1    =  1.0 )
      PARAMETER ( MAXC   =    5 )
      PARAMETER ( MAXX   =    5 )
      PARAMETER ( NKOM   =   15 )
      PARAMETER ( VLARGE = 1.E+6)
      PARAMETER ( VSMALL = 1.E-6)
      PARAMETER ( ZERO   =  0.0 )
C
      DIMENSION C      (MAXC  )
      DIMENSION CSTAR  (MAXX  ,MAXC  )
      DIMENSION G      (MAXX  )
C
      COMMON   /CVMAC / D      (MAXC  )
     1   ,              VLAM   (MAXC  )
     1   ,              VMU    (MAXC  )
      COMMON   /CVMAX1/ A      (MAXX  ,MAXX  )
     1   ,              BDL    (MAXX  )
     1   ,              BDU    (MAXX  )
     1   ,              CINVRS (MAXX  ,MAXC  )
     1   ,              CN     (MAXX  ,MAXC  )
     1   ,              DELTA  (MAXX  )
     1   ,              GM     (MAXX  )
      COMMON   /CVMAX2/ GHAT   (MAXX  )
     1   ,              H      (MAXX  ,MAXX  )
      COMMON   /KOMVMA/ KOMVMA (NKOM  )
      COMMON   /KVMAX / IBASIS (MAXX  )
C
      EQUIVALENCE ( CSTAR      ,  CINVRS ( 1,1) )
      EQUIVALENCE ( NC         ,  KOMVMA (   5) )
      EQUIVALENCE ( NVAR       ,  KOMVMA (   6) )
      EQUIVALENCE ( NVAR2      ,  KOMVMA (   7) )
C
C***********************************************************************
C
C...INITIALIZE NECESSARY VALUES
C
      IF (ISTAT .LE. -1) THEN
        CALL STINIT (ISTAT,M,MEQ,N,VLARGE)
      ENDIF
C
      DO 10 I=1,N
        GM(I) = -G(I)
   10 CONTINUE
C
      IF (M .GE. 1) THEN
        DO 20 I=1,M
          IF (I .LE. MEQ .OR. C(I) .LT. ZERO) THEN
            D(I) = ZERO
            CN(NVAR,I) = C(I)
          ELSE
            D(I) = -C(I)
            CN(NVAR,I) = ZERO
          ENDIF
C
          VLAM(I) = ZERO
   20   CONTINUE
      ENDIF
C
      BDU(NVAR) = FP1
      IFEAS = -1
   30 CONTINUE
C
C...THE QUADRATIC PROGRAMMING ROUTINE IS CALLED
C
      CALL QUAD (ISTAT,MEQ,IPRINT)
C
C...DETERMINE IF THE REQUIRED FEASIBILITY CONDITIONS HOLD
C
      IF (DELTA(NVAR) .LE. VSMALL) THEN
        ISTAT = 8
        CALL VMAMSG(ISTAT,0,ZERO)
        GO TO 999
      ENDIF
C
      DO 40 I=1,NC
        IF (IBASIS(I) .LT. NVAR2) THEN
          ISTAT = 9
          CALL VMAMSG(ISTAT,0,VLARGE)
          GO TO 999
        ELSEIF (IBASIS(I) .EQ. NVAR2) THEN
          IFEAS = 1
        ENDIF
   40 CONTINUE
C
      IF (IFEAS .EQ. 0) THEN
        ISTAT = 8
        CALL VMAMSG (ISTAT,0,ZERO)
        GO TO 999
      ELSEIF (IFEAS .LT. 0) THEN
        BDU(NVAR) = FEASP*DELTA(NVAR)
        IFEAS = 0
        GO TO 30
      ENDIF
C
C...CALCULATE THE LAGRANGE MULTIPLIERS
C
      DO 60 I=1,NC
        K = IBASIS(I) - NVAR2
        IF (K .GT. 0) THEN
          DO 50 J=1,N
            VLAM(K) = VLAM(K) + CSTAR(I,J)*GHAT(J)
   50     CONTINUE
        ENDIF
   60 CONTINUE
C
  999 CONTINUE
C
C***********************************************************************
C
      RETURN
      END
